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EXECUTIVE  SUMMARY 


The  objective  of  this  research  is  to  develop  a  groundwater  quality  modeling 
advisory  system  for  use  in  investigating  possible  remediation  activities  for  the  cleanup 
of  contamination  from  hazardous  substances,  pollutants  and  contaminants  at  Air  Force 
sites.  In  addition,  this  research  explores  the  use  of  optimization  methods  for 
determining  optimal  remediation  for  implementation  at  a  specific  site.  A  1987 
Executive  Order  authorized  the  Secretary  of  Defense  to  implement  the  Department  of 
Defense  Installation  Restoration  Program  (IRP).  The  objectives  of  this  program  include 
.  the  identification,  investigation,  research  and  development,  and  cleanup  of  sites 

contaminated  with  hazardous  substances  from  past  and  present  activities.  The  Air 
Force  has  established  its  own  in-house  management  and  technical  expertise  for 
«  implementing,  monitoring  and  managing  activities  within  the  IRP.  The  remedial  action 

process  consists  of  four  discrete  processes.  These  include:  (1)  Preliminary 
Assessment  and  Site  Inspection,  (2)  Remedial  Investigation  and  Feasibility  Study,  (3) 
Remedial  Design  and  Remediation  Action  and  (4)  Site  Closeout.  The  focus  of  the 
research  summarized  in  this  report  impacts  the  implementation  of  the  Remedial 
Investigation  and  Remedial  Action  phases  of  the  remedial  action  process. 

Over  the  past  several  decades,  many  different  models  for  contaminant 
transport  in  porous  media,  under  varying  conditions  and  assumptions,  have  been 
proposed  and  tested.  These  range  from  very  simple  models  based  on  one-dimensional 
analytical  solutions,  which  assume  a  completely  homogeneous  and  isotropic  medium, 
to  very  complex  models  based  on  three-dimensional  numerical  solutions  which  allow 
for  complete  specification  of  the  aquifer  and  contaminant  characteristics  throughout 
a  three-dimensional  grid.  All  contaminant  transport  models,  regardless  of  the 
complexity  of  the  solution  method,  require  certain  assumptions  regarding  the  nature 
of  the  transport  processes,  and,  therefore,  can  only  approximate  the  actual  spread  of 
contaminants  from  a  given  site  and  the  associated  risks  from  human  exposure  to 
contaminated  groundwater. 

This  situation  presents  a  familiar,  yet  difficult  problem  to  the  analyst  and  the 
decision-makers.  Sufficient  data  on  the  hydrogeology  are  rarely,  if  ever,  available  to 
apply  the  most  complex,  three-dimensional  contaminant  transport  models  to  a 
proposed  or  monitored  site.  The  analyst  must  choose  a  transport  model  based  on  a 
«  tradeoff  between  the  presumed  greater  accuracy  of  complex  models  and  the  less 

onerous  data  requirements  and  easier  application  of  simpler  models.  The  topic  of 
choosing  an  appropriate  model  is  one  of  the  important  aspects  of  the  advisory  system 
under  development  for  the  Air  Force,  and  specific  algorithms  have  been  developed  to 
assist  the  user  with  this  task. 

Even  with  the  choice  of  an  appropriate  transport  model,  considerable 
uncertainty  is  likely  to  be  present  in  the  analysis  of  contamination  risk.  Application 
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of  groundwater  transport  models  requires  estimation  of  parameters  which  are  both 
difficult  to  measure  and  spatially  variable,  such  as  hydraulic  conductivity  and 
dispersivity.  There  is  often  good  reason  to  doubt  the  accuracy  of  the  input  data.  For 
instance,  if  an  analytical  model  requires  the  spatial  average  of  the  hydraulic 
conductivity  throughout  the  local  area  of  the  aquifer,  and  the  available  data  consist 
of  only  one  or  two  slug  tests,  plus  perhaps  an  expert  opinion,  there  is  good  reason  to 
doubt  that  the  reported  best  estimate  of  the  parameter  accurately  reflects  the  true 
mean  value.  Simply  running  the  model  in  a  deterministic  mode  using  the  best 
estimates  of  the  parameters  may  not  provide  sufficient  information  for  a  decision, 
because  the  uncertainty  in  the  analysis  has  not  been  taken  into  account.  For  instance, 
if  a  deterministic  application  suggests  no  risk  of  contamination,  no  information  is 
provided  as  to  the  certainty  of  this  conclusion. 

The  recommended  alternative  is  to  explicitly  consider  the  uncertainty  in  the 
analysis,  through  the  use  of  Monte  Carlo  analysis.  Uncertainty  enters  the  modeling 
process  in  three  ways:  (1)  through  natural  parameter  variability;  (2)  through 
measurement  error,  which  also  introduces  uncertainty  in  parameter  estimation;  and  (3) 
through  model  error,  representing  uncertainty  introduced  by  the  degree  to  which  the 
simplifying  assumptions  used  to  develop  a  model  fail  to  accurately  represent  the  actual 
physical  processes  at  the  site  in  question.  The  first  two  of  these  sources  of 
uncertainty  can  be  analyzed  separately.  However,  the  data  are  often  insufficient:  in 
such  cases,  the  natural  and  measurement  uncertainty  may  be  combined  into  one 
source  of  uncertainty  for  the  Monte  Carlo  analysis,  through  the  specification  of  the 
distribution  of  the  parameter  value. 

The  third  source  of  uncertainty  in  the  analysis  is  due  to  the  degree  to  which  the 
transport  model  applied  may  misrepresent  actual  processes  at  the  site.  Examples  of 
this  source  of  uncertainty  include  the  sorption  of  contaminants  to  soil  surfaces  and 
degradation  rate  coefficients.  This  source  of  uncertainty  is  very  difficult  to  quantify, 
and  indeed  may  be  impossible  to  quantify  for  specific  sites,  unless  extensive  sampling 
and  monitoring  data  are  available. 

A  Monte  Carlo  analysis  requires  that  distributions  be  specified  for  the  underlying 
parameters  having  the  greatest  impact  on  contaminant  transport.  Specification  of  a 
parameter  distribution  consists  of  two  steps:  (1)  choice  of  a  distributional  form,  and 
(2)  specification  of  the  descriptive  parameters  of  that  distribution.  On  the  first  issue, 
the  choice  of  distributional  form,  the  system  does  of  necessity  provide  some 
limitations.  That  is,  for  models  which  are  expected  to  be  used  in  cases  for  which  the 
impacted  aquifer  is  at  least  moderately  well-characterized,  certain  parameter 
distributions  are  constrained  to  follow  specific  forms,  which  are  generally  well 
accepted  in  the  literature.  For  instance,  in  some  of  the  models  the  mean  hydraulic 
conductivity  must  be  specified  by  a  log-normal  distribution.  However,  even  in  these 
cases,  a  choice  is  present  in  the  parameterization,  as  the  mean  hydraulic  conductivity 
may  be  directly  specified  from  the  log-normal,  or  generated  from  underlying  parameter 


distributions.  In  general,  where  the  parameters  are  at  least  moderately  well  known  the 
choice  of  a  distributional  form  should  not  have  a  major  impact  on  the  results.  In  its 
present  form,  the  Advisory  System  incorporates  the  framework  for  Monte  Carlo 
analysis,  but  additional  research  is  needed  to  develop  the  parameter  distributions  and 
values  for  site-specific  hydrogeologic  conditions. 

in  addition  to  aiding  in  the  choosing  of  an  appropriate  mathematical  model  for 
a  specific  site,  the  Advisory  System  is  being  modified  to  determine  efficient  or 
optimal  remediation  strategies.  The  optimization  routine  evaluates  tradeoffs  between 
the  long-term  cost  of  remediation  and  the  probability  the  remediation  strategy  will  fail. 

A  chance-constrained  optimization  model  is  being  developed  to  determine  the 
most  efficient  groundwater  remediation  strategies.  The  optimization  model  is 
multiobjective  and  driven  by  probabilistic  measures  of  contaminant  concentration  in 
the  groundwater  surrounding  the  hazardous  waste  site.  The  chance-constrained 
model  is  used  to  determine  the  tradeoffs  that  exist  between  short-term  and  long-term 
remediation  costs  and  the  probability  that  the  remediation  strategy  will  fail. 

The  development  of  an  efficient,  effective  and  reliable  remediation  strategy 
requires  a  clear  understanding  of  the  site  characteristics  and  the  remediation  actions 
implemented.  In  addition,  the  optimal  remediation  strategy  must  consider  tradeoffs 
between  the  remediation  cost  and  the  reliability  of  the  remediation  strategy.  By 
investigating  these  tradeoffs,  the  decision  maker  can  more  accurately  assess 
remediation  needs,  feasible  remediation  strategies  and  remediation  strategy 
effectiveness. 

Long-term  remediation  costs  depend  on  specific  remediation  considerations  and 
actions.  Examples  of  possible  remediation  strategies  include  pulse  pumping  and 
treatment,  and  continuous  pumping  and  treatment.  Potential  cost  savings  are  realized 
by  varying  the  long-term  remediation  action.  The  reliability  of  the  long-term 
remediation  strategy  represents  the  likelihood  that  contaminant  concentrations  within 
the  groundwater  exceeds  specified  maximums  and  are  modeled  as  constraints.  These 
two  conflicting  goals  or  objectives  are  weighed  against  one  another,  using  a  chance- 
constrained  optimization  model  in  which  the  physical  constraints  are  originally 
expressed  as  probabilistic  statements. 

Using  this  methodology,  optimal  groundwater  remediation  strategies  are 
determined  by  minimizing  the  long-term  and  short-term  costs  associated  with  the  site 
remediation.  In  addition,  the  optimal  remediation  strategies  are  conditioned  on  the 
probability  that  the  contaminant  concentration  at  any  time  does  not  exceed 
pre-specified  maxima.  The  actual  concentrations  at  any  specified  coordinates  are 
calculated  by  solving  the  governing  differential  equations  for  groundwater  contaminant 
flow  in  which  key  site  characteristics  are  expressed  as  random  variables.  The 
resulting  optimization  model  is  solved  using  a  second  moment  formulation  combined 
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with  Monte  Carlo  simulation. 


The  final  product  of  this  project  is  a  computer-based  Air  Force  Installation 
Restoration  Advisory  System  Workstation  for  contaminant  modeling  and  decision 
making.  When  completed,  the  Advisory  System  will  be  fully  documented  and 
compatible  with  the  DOS  and  UNIX  operating  systems.  This  software  can  be  used  as 
an  aid  to  technical  project  managers  within  the  U.S.  Air  Force  Installation  Restoration 
Program  in  developing  and  evaluating  possible  remediation  alternatives  and  managing 
ongoing  remediation  activities. 
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SECTION  I 


INTRODUCTION 


A.  OBJECTIVES 

Air  Force  needs  for  contaminant  transport  mathematicai  modeling  and  decision¬ 
making,  in  terms  of  the  predictive  requirements  of  the  Installation  Restoration  Program 
(IRP),  may  at  least  be  partially  addressed  by  development  of  an  interactive,  user- 
friendly  computer-based  engineering  workstation.  Background  information  on 
Department  of  Defense  environmental  restoration  efforts  and,  specifically,  the  Air 
Force  IRP  is  presented  in  the  next  section.  The  principal  element  of  the  workstation 
is  an  advisory  system  incorocrating  basic  software  to:  define  the  magnitude,  extent, 
direction,  and  rate  of  movement  of  identified  contaminants;  identify  significant  public 
health  and  environmental  hazards  of  migrating  pollutants;  recommend  candidate 
remedial  actions;  maintain  databases  of  model  parameters,  and  accomplish  other 
supporting  tasks.  The  principal  function  of  a  workstation  is  to  provide  optimal  and 
efficient  support  to  its  user  regarding  the  tasks  determined  for  the  user/workstation 
entity.  Generally,  this  function  can  be  divided  into  a  number  of  subfunctions  which 
are  determined  by  analyzing  the  tasks  performed  by  the  intended  users  and  the 
hardware/software  environments  available.  Important  elements  in  this  analysis  are 
determining  the  amount  and  type  of  data  and  establishing  the  level  of  synthesis 
required  to  adequately  perform  the  required  tasks.  Furthermore,  the  various  levels  of 
expertise  of  potential  users  must  be  determined  and  accommodated  for  in  the 
operating  system  to  provide  users  with  adequate  assistance. 


B.  BACKGROUND 

The  legal  mandate  for  the  Air  Force  (AF)  Installation  Restoration  Program  (IRP) 
is  the  Comprehensive  Environmental  Response,  Compensation  and  Liability  Act  of 
1 980  (CERCLA,  known  as  the  Superfund  Act)  and  the  Superfund  Amendments  and 
Reauthorization  Act  of  1986  (SARA).  Section  21 1  of  SAFIA  deals  with  the  Defense 
Environmental  Restoration  Program  (DERP),  of  which  the  IRP  is  the  primary 
subcomponent  (Reference  1).  A  1987  Executive  Order  provided  authority  to  the 
Secretary  of  Defense  to  implement  the  Department  of  Defense  (DOD)  Environmental 
Restoration  Program  within  the  overall  framework  of  CERCLA  and  SARA.  The 
objectives  of  the  IRP  include  "the  identification,  investigation,  research  and 
development,  and  cleanup  of  contamination  from  hazardous  substances,  pollutants, 
and  contaminants.”  The  program  is  focused  on  cleanup  of  detected  contamination 
from  past  activities,  but  as  noted  includes  research  as  well  as  development  and 
demonstration  of  innovative  and  cost-effective  cleanup  technologies.  IRP  activities  are 
managed  centrally  in  the  Office  of  the  Secretary  of  Defense  and  are  carried  out  by  the 
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Military  Services  and  Defense  Agencies.  Under  this  agreement,  the  U.S.  Air  Force 
retains  the  authority  and  initiative  for  cleanup  activities  at  its  own  installations. 


The  Air  Force  has  established  its  own  in-house  management  and  technical 
expertise  for  implementing  the  IRP,  following  a  decentralized  approach  which  places 
emphasis  and  authority  with  the  Major  Air  Commands  (MAJCOMs)  and,  in  turn,  with 
the  individual  installations  under  their  jurisdiction  (Reference  2).  Several  service 
organizations  support  the  implementation  of  the  Air  Force  IRP:  the  Air  Force  Civil 
Engineering  Support  Agency  (AFCESA,  HQ  at  Tyndall  AFB,  Florida),  Armstrong 
Laboratory  Environics  Directorate,  Tyndall  AFB,  Florida,  the  Center  for  Environmental 
Excellence  (AFCEE,  Brooks  AFB,  San  Antonio,  Texas),  and  the  AF  Regional  Civil 
Engineer  offices.  Additional  support  is  provided  by  the  Air  Force  Material  Command 
(AFMC),  which  is  responsible  for  the  advancement  and  effective  management  of  the 
Air  Force  scientific  and  technical  resources.  An  Air  Force  Installation  Restoration 
Management  (AFIRM)  Committee  has  also  been  organized  to  support  the  MAJCOMs 
and  review  remedial  action  plans  for  complex  problems. 

The  remedial  action  process  is  a  progression  of  steps  designed  to  fully  analyze 
and  address  site  problems,  grouped  functionally  by  stages,  as  follows: 

1.  Preliminary  Assessment/Site  Inspection  (PA/SI)  Stage, 

2.  Remedial  Investigation/Feasibility  Study  (RI/FS)  Stage, 

3.  Remedial  Design/Remedial  Action  (RD/RA)  Stage, 

4.  Site  Closeout  (SC)  Stage. 

Figure  1  illustrates  these  four  stages  and  14  steps  of  the  remedial  action 
process.  The  opportunity  for  application  of  contaminant  transport  models  arises 
primarily  in  the  second  (investigation)  and  third  (cleanup)  stages.  However, 
mathematical  models  may  be  used  in  the  first  stage  in  the  case  of  unknown 
subsurface  sources  of  contamination:  the  most  likely  location  of  the  source  could  be 
calculated  from  known  field  measurements  of  the  edge  of  the  plume  —  as  part  of  the 
discovery  and  preliminary  assessment  steps. 

in  the  second  stage,  mathematical  models  may  be  applied  to: 

•  Estimate  the  rate  and  extent  of  contamination  migration  from  several 
sources  (surface  and  subsurface); 

•  Simulate  current  and  future  scenarios  of  contamination  and  potential 
impacts  at  all  locations  of  interest; 

•  Evaluate  the  likely  effectiveness  of  proposed  alternatives  for  remediating 
the  impacts  of  released  contaminants; 
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•  Perform  risk  analysis,  accounting  for  uncertainty  in  predictions,  to  assist 
in  selection  of  the  best  remedial  strategy. 


Figure  1.  IRP  Remedial  Action  Process 


In  the  third  stage  (cleanup),  models  are  useful  in  designing  the  remedial 
strategy:  the  optimal  strategy  should  be  cost-effective.  Models  do  not  reduce  the 
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need  for  good  quality  site*specific  data:  they  help  determine  data  needs,  make  better 
use  of  available  data,  and  refine  the  data  collection  (monitoring)  process  to  insure 
compliance  with  cleanup  goals. 

The  Air  Force  Center  for  Environmental  Excellence  (AFCEE)  at  Brooks  AFB  (San 
Antonio,  Texas)  operates  the  technical  information  management  system  (IRPIMS)  for 
Air  Force  IRP  sites.  It  is  one  of  the  contract  support  centers  for  investigative  studies. 
It  can  provide  technical  consultation,  field  monitoring,  sample  analysis  support,  and 
has  developed  programs  on  site  ranking  and  Quality  Assurance/Quality  Control 
(QA/QC). 


C.  SCOPE 

Seven  groundwater  models  and  one  surface  water  quality  model  are  maintained 
within  IRPIMS.  However,  these  models  need  updating,  do  not  have  the  capability  of 
addressing  uncertainty  in  model  prediction,  are  not  integrated  within  a  decision- 
analysis  framework,  and  do  not  address  many  of  the  complex  flow  and  mass  transport 
phenomena  exhibited  by  solvents  in  groundwater. 

Advisory  systems  provide  a  systematic  framework  for  the  determination  of 
optimal  choices  independent  of  (but  interacting  with)  the  transport  models  employed 
for  concentration  calculations.  The  conceptual  structure  of  such  systems  incorporates 
several  premises: 

1 .  The  best  decisions  are  made  by  technically  informed  individuals  which 
may  be  aided  by  computer-based  data  handling  and  analysis, 

2.  Synthesis  of  technical,  socio-economic  and  political  aspects  may  be 
necessary; 

3.  Problem  visualization  through  high-resolution  graphics  facilitates  insight 
and  problem  comprehension; 

4.  Guidance  and  problem  representation  assistance,  not  dictated  decisions, 
provide  optimal  decision  support; 

5.  The  central  role  of  the  human  decision  maker  requires  that  decision 
support  tools  be  designed  as  man-machine  systems. 

A  workstation  may  include  both  surface  and  groundwater  models,  aimed  at 
developing  alternative  remediation  strategies  for  polluted  surface  and  groundwater 
systems,  and  at  designing  the  technical  details  of  a  preferred  remedial  action.  It 
should  be  designed  to  perform: 
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1 .  Data  management  and  analysis 

a.  Acquisition  of  data  directly  from  the  field,  from  maps, 
and  from  other  databases: 

b.  Data  reduction; 

c.  Data  analysis: 

d.  Data  storage  and  retrieval;  and 

e.  Graphical  presentation  of  data. 

2.  Site  characterization 

3.  Source  identification 

4.  Plume  delineation 

5.  Contaminant  transport  analysis 

6.  Risk  analysis 

7.  Evaluation  and  optimization  of  potential  remedial  action  alternatives, 
compliance  monitoring,  sampling  strategies 

Its  components  should  include: 

1 .  Hardware 

a.  Graphics  capability 

b.  Peripherals  (e.g.,  printer,  mouse) 

c.  Communication  links 

d.  Storage  devices 

2.  Software 

a.  Data  storage,  management,  analysis 

b.  Simulation  models 

c.  Hydrogeologic,  hydrogeochemical  analysis 

d.  Shell:  user-friendly  interface 

I.  help  screens 

ii.  advisory  system  interface 

iii.  database  of  model  parameters 

iv.  linkage  to  IRPIMS,  graphics 

In  designing  the  workstation,  flexible  architecture  is  necessary  for  efficient 
updating,  maintenance,  and  expansion  of  hardware  and  software,  in  addition,  an 
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operational  support  structure  needs  to  be  implemented  for  the  system  maintenance 
and  to  provide  user-application  support.  Finally,  as  an  integral  part  of  the 
organizational  workstation  environment,  a  continuing  technology  transfer  program 
should  be  developed  to  include  general  introduction  courses,  various  levels  of  on-site 
hands-on  training,  and  roving  experts  visiting  the  different  workstation  locations  on 
a  regular  basis. 
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SECTION  II 


METHODOLOGY 


A.  ROLE  OF  UNCERTAINTY  IN  GROUNDWATER  RISK  ANALYSIS 

Over  the  past  several  decades  many  different  mathematical  models  for  contai 
inant  transport  in  porous  media  have  been  proposed  and  tested,  under  varying 
conditions  and  assumptions.  These  range  from  models  based  on  closed-form 
analytical  solutions  to  one-,  two-  and  three-dimensional  versions  of  the  governing 
differential  equations  (assuming  completely  homogeneous  and  isotropic  media),  to 
highly  complex  numerical  solutions  of  equations  which  allow  for  complete 
specification  of  both  the  aquifer  and  the  contaminant  characteristics  throughout  a 
three-dimensional  grid.  Ail  of  these  contaminant  transport  models,  regardless  of  the 
complexity  of  the  solution  method,  require  certain  assumptions  regarding  the  nature 
of  the  transport  processes  and  the  physical  system  abstracted,  and  so  can  provide 
only  an  approximation  of  the  actual  spread  of  the  contaminant(s)  from  a  given  site  and 
the  associated  risk. 

This  situation  presents  a  familiar,  yet  difficult  problem  to  the  analyst  and  the 
decision-makers.  Sufficient  data  on  the  hydrogeology  are  rarely  available  to  apply  the 
most  complex,  three-dimensional  flow  and  contaminant  mass  transport  models  to 
even  well  monitored  sites.  The  analyst  must,  whether  explicitly  or  implicitly,  choose 
a  transport  model  based  on  a  trade-off  between  the  presumed  greater  accuracy  of 
complex  models  and  the  less  onerous  data  requirements  and  easier  application  of 
simpler  models.  The  topic  of  choosing  an  appropriate  model  is  one  of  the  important 
aspects  of  the  advisory  system  under  development  for  the  Air  Force. 

Even  with  the  choice  of  an  appropriate  code-verified  transport  model,  consid¬ 
erable  uncertainty  is  likely  to  be  present  in  the  analysis  of  contamination  risk.  In 
groundwater  transport  models  (which  require  estimation  of  parameters  which  are 
difficult  to  measure  and  spatially  variable,  such  as  hydraulic  conductivity  and 
dispersivity),  there  is  often  good  reason  to  doubt  the  accuracy  of  the  input  data.  For 
instance,  if  an  analytical  model  requires  the  spatial  average  of  the  hydraulic 
conductivity  throughout  the  local  area  of  the  aquifer,  and  the  available  data  consist 
of  only  one  or  two  slug  tests,  plus  perhaps  an  expert  opinion,  there  is  good  reason  to 
doubt  that  the  reported  best  estimate  of  the  parameter  accurately  reflects  the  true 
mean  value.  Simply  running  the  model  in  a  deterministic  mode  using  the  best 
estimates  of  the  parameters  will  not  then  provide  sufficient  information  for  a  decision, 
because  the  uncertainty  in  the  analysis  has  not  been  taken  into  account.  For  instance, 
if  a  deterministic  application  suggests  no  risk  of  contamination,  no  information  is 
provided  as  to  the  certainty  of  this  conclusion. 
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The  recommended  alternative  is  to  explicitly  consider  the  uncertainty  that  is 
present  in  the  analysis,  through  the  use  of  stochastic  methods  such  as  Monte  Carlo 
analysis.  Uncertainty  enters  the  modeling  process  in  three  ways:  (1 )  through  natural 
parameter  variability;  (2)  through  measurement  error,  which  also  introduces 
uncertainty  in  parameter  estimation;  and  (3)  through  model  error,  representing 
uncertainty  introduced  by  the  degree  to  which  the  simplifying  assumptions  used  to 
develop  a  model  fail  to  accurately  represent  the  actual  physical  processes  at  the  site 
in  question.  The  first  two  of  these  sources  of  uncertainty  can  be  analyzed 
separately.  However,  the  data  are  often  insufficient:  in  such  cases  the  natural  and 
measurement  uncertainty  may  be  combined  into  one  source  of  uncertainty  for  the 
Monte  Carlo  analysis,  through  the  specification  of  the  distribution  of  the  parameter 
value. 


The  third  source  of  uncertain;.  the  analysis  is  due  to  the  degree  to  which  the 
transport  model  applied  may  misrepresent  actual  processes  at  the  site.  This  source 
of  uncertainty  is  unfortunately  very  difficult  to  quantify,  and  may  indeed  be  impossible 
to  do  so  for  specific  sites  unless  extensive  sampling  and  monitoring  data  are 
available. 

To  conduct  a  Monte  Carlo  analysis,  it  is  required  that  distributions  be  specified 
for  the  underlying  parameters.  Specification  of  a  parameter  distribution  consists  of 
two  steps:  choice  of  a  distributional  form,  and  specification  of  the  hyperparameters 
of  that  distribution.  On  the  first  issue,  the  choice  of  distributional  form,  the  system 
does  of  necessity  provide  some  limitations.  That  is,  for  models  which  are  expected 
to  be  used  in  cases  for  which  the  impacted  aquifer  is  at  least  moderately  well 
characterized,  certain  of  the  parameter  distributions  are  constrained  to  follow  specific 
forms,  which  are  generally  well  accepted  in  the  literature.  For  instance,  in  some  of 
the  models  the  mean  hydraulic  conductivity  must  be  specified  by  a  log-normal 
distribution.  However,  even  in  these  cases  a  choice  is  present  in  the  parameterization, 
as  the  mean  hydraulic  conductivity  may  be  directly  specified  from  the  log-normal,  or 
generated  from  underlying  parameter  distributions.  In  general,  where  the  parameters 
are  at  least  moderately  well  known  the  choice  of  a  distributional  form  should  not  have 
a  major  impact  on  the  results.  In  its  present  form,  the  workstation  Advisory  System 
incorporates  the  framework  for  Monte  Carlo  analysis  for  several  models,  but  additional 
research  is  needed  to  develop  the  required  parameter  distributions  and  values  for  site- 
specific  hydrogeologic  conditions. 

Once  a  specific  transport  model  is  selected,  an  estimate  of  the  unconditional 
distribution  of  the  contaminant  concentration  is  needed  to  assess  the  risk  associated 
with  a  site.  The  resulting  unconditional  contaminant  concentration  distribution  is  then 
utilized  to  assess  the  risk  associated  with  a  specific  site  using  Bayes'  theorem.  For 
example,  if  variabilities  are  considered  in  the  time  to  the  waste  container  failure,  the 
actual  leachate  release  concentration  and  the  spatial  variation  of  the  hydraulic 
conductivity,  the  unconditional  distribution  of  the  contaminant  concentration  for  any 
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location  and  time  may  be  expressed  as  (References  3,  4,  5,  6  ): 


PlCixyt)]  -  /////  [  PiCixyt)  \Cj,,  t^,K)  p{Cj\tf)  • 
p{t/|p}  p{p^  p{jC|Ujf}  p{ujj|vjf}  P  (vjf)  ]  dC^  dtf  dp  dv*  du* 


where  ff  »  time  to  containment  failure  (or  initiation  of  contamination),  Q  »  leachate 
release  concentration,  K  »  log  hydraulic  conductivity  field,  p  =  parameter  of  the 
distribution,  u  »  mean  of  the  distribution,  and  v  =  variance  of  the  distribution.  In  the 
deterministic  case,  the  contaminant  concentration  may  be  determined  conditioned  on 
the  specific  combination  of  parameter  values  chosen.  However,  in  reality,  each  of  the 
model  parameter  values  may  be  considered  to  be  a  random  variable  with  a  specified 
distribution  or  a  mean,  variance  and  an  assumed  distribution.  In  this  case,  the 
contaminant  concentration  must  be  evaluated  as  a  function  of  each  random  variable. 
Solution  of  the  above  integral  is  intractable,  but  an  alternative  is  combining  probability 
distributions  of  input  parameters  with  a  deterministic  model.  Such  Monte  Carlo 
simulation  experiments  were  conducted  by  Medina  et  al.  (References  5,  6)  with  a  two 
dimensional  numerical  solute  transport  model.  The  effect  on  the  probability 
distribution  of  the  contaminant  concentration  (at  an  observation  point  in  the  flow  field) 
due  to  variance  in  the  input  parameters,  is  illustrated  in  Figure  2.  There  is  only  one 
concentration  associated  with  a  probability  of  1.0  for  the  deterministic  solution, 
whereas  a  spread  in  the  probability  distributions  represents  the  uncertainty  in  the 
predictions  where  variability  in  the  input  parameters  is  allowed.  Yet,  the  uncertainty 
in  the  predictions  is  over  a  much  broader  range  about  the  deterministic  solution  for 
variance  in  hydraulic  conductivity  than  for  variance  in  time  to  failure. 

A  major  problem  in  determining  the  risks  of  any  site  are  the  uncertainties 
associated  with  model  parameters  such  as  the  leachate  release  concentration,  the 
hydraulic  conductivity  field  and  the  time  to  containment  failure.  In  most  cases  only 
mean  and  variance  of  the  distributions  of  the  individual  parameters  are  known.  To 
incorporate  model  parameter  uncertainties,  approximate  solution  techniques  may  be 
used  (e.g.,  Monte-Carlo  simulation).  An  estimate  of  the  cumulative  distribution 
function  for  the  unconditional  contaminant  concentration  is  then  developed.  The 
unconditional  distribution  can  then  be  used  to  assess  the  risks  due  to  the 
contamination  at  any  site.  Using  Bayes'  theorem,  a  more  realistic  distribution  for  the 
contaminant  concentrations  may  be  developed  conditioned  on  additional  input  data. 
Additional  site  information  concerning  any  of  the  model  parameters  would  result  in  a 
reduction  in  the  uncertainty  of  that  parameter,  which  could  result  in  a  reduction  in  the 
error  associated  with  predicting  the  groundwater  contaminant  concentration. 

To  begin,  let  the  uncertainty  associated  with  the  random  variable  A  be  defined 
by  the  prior  distribution  P(A).  The  effect  of  additional  information  collected  through 
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Figure  2.  Uncertainty  in  Subsurface  Contaainant  Modeling 
Due  to  Variance  in  Input  Parameters 
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field  observation  and  characterized  by  the  conditional  probability  distribution,  PdtSA), 
can  now  be  defined  using  Bayes'  theorem: 


pwii,,  ._p<£i!fLzw>_ 

J  P(Ij^A)  P(A)  dA 


where  P(A\IJ  is  the  updated  distribution  of  the  random  variable  A  and  is  referred  to 
as  the  posterior  distribution.  The  posterior  distribution  reflects  the  updated 
distribution  of  the  random  variable  A. 

Bayes'  theorem  provides  an  ideal  tool  for  continually  assessing  the  errors  and 
concentrations  associated  with  groundwater  contaminant  transport.  A  probabilistic 
assessment  of  the  groundwater  contaminant  concentration  and  its  corresponding  error 
can  then  be  used  to  address  the  risks  associated  with  a  specific  site.  Regulatory 
actions  or  remedial  decisions  based  on  this  approach  can  be  significantly  different  and 
more  realistic  from  those  based  on  a  deterministic  estimate  of  groundwater 
contaminant  concentrations  when  a  large  amount  of  variance  is  present  in  the 
unconditional  distribution  of  the  contaminant  concentration  (Reference  6). 

This  approach  has  the  distinct  advantage  of  being  able  to  more  accurately 
assess  the  groundwater  contaminant  concentrations  at  a  specific  site.  In  addition, 
this  methodology  allows  for  the  immediate  incorporation  of  new  observational  data 
for  the  site  in  an  effort  to  further  refine  the  assessment.  This  decision-making 
approach  could  be  adapted  for  Air  Force  IRP  needs. 


B.  OVERVIEW  OF  THE  ADVISORY  SYSTEM 

A  flow  chart  illustrating  the  design  of  the  workstation  Advisory  System  is 
presented  in  Figure  3.  The  user/anaiyst  interacts  with  a  module  that  controls  the  flow 
between  the  various  elements  of  the  system.  For  example,  to  the  left  of  the  system 
manager  are  modules  that  access  stored  data  (site-specific  data,  regional  data,  data 
on  model  input  parameters)  and  preliminary  screening  modules  (to  rank  the  severity 
of  contamination  at  the  site  under  investigation).  To  the  right  of  the  manager  module 
is  a  transport  model  selection  module,  named  the  CHOICE  algorithm,  discussed  in 
much  greater  detail  in  a  later  sub-section.  It  essentially  aids  the  inexperienced  user 
in  selection  of  the  solute  transport  model  most  appropriate  for  the  site  hydrogeology 
and  method  of  waste  disposal.  After  the  appropriate  selection  is  made  a  plume  is 
predicted  and  a  cumulative  probability  distribution  of  contaminant  concentration  is 
derived  at  any  desired  point  in  the  flow  field.  The  amount  of  variance  in  the  prediction 
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Figur*  3.  Workstation  Advisory  System  Flow  Chart 


indicates  the  degree  of  uncertainty,  which  can  be  reduced  by  additional  field  sampling. 
The  next  step  is  optimizing  the  remediation  process,  providing  a  framework  for 
evaluating  remediation  alternatives  and  implementing  a  solution  at  minimal  cost  and 
environmental  risk.  An  algorithm  to  select  remediation  alternatives  is  currently  under 
development.  Details  of  the  optimization  process  are  presented  below.  If  further 
relevant  field  data  is  available,  the  cycle  of  transport  modeling  begins  again,  to 
possibly  reduce  the  variance  in  the  predictions. 


C.  OPTIMIZATION  OF  REMEDIATION 

The  chance  constrained  optimization  model  is  multi-objective  in  nature  and 
driven  by  probabilistic  measures  of  contaminant  concentration  in  the  groundwater 
surrounding  the  hazardous  w^ste  s’te.  Chance  constraints  are  used  to  determine  the 
trade-offs  that  exist  between  short-term  and  long-term  remediation  costs  and  the 
probability  that  the  remediation  strategy  will  fail.  It  is  important  to  emphasize  that  the 
whole  process  depends  on  the  transport  models  (described  in  detail  below),  which  are 
executed  in  probabilistic  mode.  The  development  of  an  efficient,  effective  and  reliable 
remediation  strategy  requires  a  clear  understanding  of  the  site  characteristics  and  the 
remediation  actions  implemented.  In  addition,  the  optimal  remediation  strategy  must 
consider  trade-offs  between  the  remediation  cost  and  the  reliability  of  the  remediation 
strategy.  By  investigating  these  trade-offs,  the  decision  maker  can  more  accurately 
assess  remediation  needs,  feasible  remediation  strategies  and  remediation  strategy 
effectiveness. 

Long-term  remediation  costs  will  depend  on  specific  remediation  considerations 
and  actions.  Examples  of  possible  remediation  strategies  include  pulse  pumping  and 
treatment,  and  continuous  pumping  and  treatment.  Potential  cost  savings  can  be 
realized  by  varying  the  long-term  remediation  action.  The  reliability  of  the  long-term 
remediation  strategy  represents  the  likelihood  that  contaminant  concentrations  within 
the  groundwater  exceeds  specified  maximums  and  can  be  modelled  as  constraints. 
These  two  conflicting  goals  or  objectives  can  be  weighed  against  one  another  using 
a  chance  constrained  optimization  model. 

Chance  constrained  modeling  is  an  optimization  method  in  which  the  physical 
constraints  are  originally  expressed  as  probabilistic  statements.  The  technique  was 
first  applied  to  heating  oil  refinement  operations  (Reference  7).  More  recently,  chance 
constrained  optimization  has  been  applied  in  water  resources  to  determine  the  optimal 
design  and  operation  of  reservoirs  conditioned  on  stochastic  streamflows  and  to 
manage  stream-aquifer  systems  (References  8,  9,  10).  Optimal  groundwater 
remediation  strategies  are  determined  by  minimizing  the  long-term  and  short-term 
costs  associated  with  the  site  remediation.  In  addition,  the  optimal  remediation 
strategies  are  conditioned  on  the  probability  that  the  contaminant  concentration  at  any 
time  does  not  exceed  pre-specified  maxima.  The  actual  concentrations  at  any 


specified  coordinates  are  calculated  by  solving  the  governing  differential  equations  for 
groundwater  contaminant  transport,  in  which  key  model  parameters  are  expressed  as 
random  variables.  The  resulting  optimization  model  can  then  be  solved  using  one  of 
several  well-established  optimization  techniques  that  include  Monte  Carlo  simulation, 
heuristic  search  techniques  and  nonlinear  optimization  techniques. 

Initially,  two  objectives  are  being  used  to  drive  the  optimization  model.  The 
first  objective  is  to  minimize  the  long-term  remediation  costs.  One  example  is  the  use 
of  waste  transport  simulation  to  locate  pumping  wells,  so  that  contamination 
concentrations  can  be  lowered  below  acceptable  health  standards  in  the  shortest 
period  of  time.  The  second  objective  will  be  to  maximize  the  reliability  of  the 
long-term  remediation  strategy  by  minimizing  the  probability  that  contaminant 
concentrations  exceed  maximum  allowable  levels  at  the  site. 

Realization  of  these  two  objectives  will  be  governed  by  chance  constraints  that 
stipulate  acceptable  levels  of  groundwater  contaminant  concentrations  at  the 
hazardous  waste  monitoring  sites.  Using  the  three  dimensional  form  of  the  solute 
transport  equation  for  uniform  groundwater  flow,  the  model  constraint  is  developed 
as  a  probability  statement  that  represents  the  likelihood  of  exceeding  specified  upper 
bounds  on  contaminant  concentrations.  The  impact  of  groundwater  treatment  and 
remediation  operations  are  directly  incorporated  into  the  governing  differential 
equations  for  contaminant  fate  and  transport  in  the  groundwater  surrounding  the  site. 
The  general  form  of  the  optimization  model  is: 

minimize  =  COSTi  C(x,y,t)) 
maximize  Z,  =  meuc  {«} 


subject  to: 

PC  C{x,y,  t)  s  MCL  ]  z  a  'i  x,y,t 
C{x,  y,  t)  z  0  x,y,t 


where  Cfx,y,t)  is  the  contaminant  concentration  at  any  point  and  time,  MCL  is  the 
specified  maximum  allowable  contaminant  concentration  and  a  is  the  desired  system 
reliability  for  any  point  in  the  subsurface  and  time  in  the  planning  horizon.  Figure  4 
presents  one  possible  algorithm  for  evaluating  the  chance  constrained  optimization 
model. 

One  major  task  of  solving  the  optimization  model  is  evaluating  the  chance 
constraint.  Due  to  the  variability  of  groundwater  parameters  and  numerical  problems 
associated  with  solving  the  governing  groundwater  equations,  the  chance  constraint 
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cannot  be  evaluated  directly.  One  alternative  is  to  use  Monte  Carlo  simulation. 
However,  typical  Monte  Carlo  simulation  requires  a  large  number  of  iterations  and  can 
become  prohibitive  with  respect  to  the  required  solution  time.  In  addition,  data 
associated  with  the  physical  parameters  is  often  lacking  for  a  specific  site  making  it 
difficult  to  fully  define  the  probability  distributions  for  important  transport  properties 
such  as  the  hydraulic  conductivity.  To  reduce  the  number  of  Monte  Carlo  iterations 
required  while  maintaining  solution  accuracy,  a  second  moment  formulation  is  used 
to  evaluate  the  chance  constraint  (Reference  11). 

To  illustrate  the  methodology,  an  example  of  preliminary  results  is  presented 
for  Operational  Unit  3  (OU3)  at  Hill  AFB,  Utah  in  Section  III.  This  example  considers 
only  pump  and  treat  activities  as  a  possible  remediation  strategy.  In  addition,  the 
hydraulic  conductivity,  retardation  by  sorption  and  porosity  are  modeled  as  random 
variables. 


D.  COMPONENT  MODELS  OF  THE  ADVISORY  SYSTEM 
1 .  CHOICE,  Algorithm  for  Model  Selection 

Ultimately,  the  management  of  any  system  means  making  decisions  aimed  at 
achieving  the  system's  goals  without  violating  specified  technical  and  nontechnical 
constraints  imposed  on  it  (Reference  1 2).  The  objective  function  is  to  minimize  costs 
and  maximize  the  effectiveness  of  remediation,  which  can  also  be  expressed  as 
minimizing  the  probability  of  failure.  This  probability  of  failure  may  be  defined  as  the 
probability  of  exceeding  a  regulatory  standard. 

The  nature  of  the  overall  modeling  process  (of  which  model  selection  is  just 
one  step)  may  be  summarized  in  five  general  steps  (Reference  13): 

•  problem  characterization  —  the  analyst  clearly  identifies  the 
exposure  assessment  study  objectives  and  constraints; 

•  site  characterization  —  the  analyst  reviews  all  available  data, 
and  possibly  develops  a  "conceptual"  model; 

•  model  selection  criteria  —  the  analyst  matches  the  objective, 
technical  and  implementation  criteria  to  available  models  and 
selects  the  most  appropriate  model(s),  in  this  case  with  the 
aid  of  the  CHOICE  algorithm; 

•  code  installation  —  in  the  case  of  a  computer  code,  the 
model(s)  should  be  properly  installed  and  tested  with  accepted 
solutions  to  standard  problems; 

•  model  application  —  the  verified  model  uses  site  data  as  input 
for  the  contaminant  assessment. 


CHOICE  is  not  a  predictive  model,  but  rather  a  screening  model.  The  algorithm 
requests  information  about  the  means  of  waste  disposal  (e.g.,  lagoons,  landfills,  rotary 
distributors,  spray  irrigation  devices,  etc.),  the  nature  of  the  aquifer,  the  perimeter  of 
compliance,  penetration,  type  of  waste  and  many  other  factors.  The  selection 
algorithm  is  part  of  an  interactive,  menu-driven  management  program  which  executes 
a  large  number  of  supporting  decision  algorithms  and  mathematical  models.  The 
mathematical  details  of  the  models  are  presented  In  the  next  sections,  and  the 
theoretical  basis  of  the  management  modules  is  presented  elsewhere  (References  3, 
5,  6,  14  and  1 5).  Criteria  for  choosing  among  transport  models  has  also  been  a  topic 
of  regulatory  interest  (References  12  and  13),  but  without  guiding  the  user  to  a 
specific  model.  The  following  contaminant  transport  models  have  presently  been 
incorporated  into  the  workstation  advisory  system: 

a.  Analytical  Models 

i.  One-Dimensional  Transport  Model  ODAST  (Reference  16). 

ii.  Two-Dimensional  Transport  Model  TDAST  (Reference  16). 

iii.  Two-Dimensional  Transport  Model  PLUM2D  (Reference  17). 

iv.  Two-Dimensional  (x,z)  Transport  Model  DUPVG 
(Reference  18). 

V.  Three-Dimensional  EPA  Monte  Carlo  Transport  Model 
EPAGW  (Reference  19). 

vi.  EPA  Monte  Carlo  Transport  Model  for  Impact  on  Surface 
Waters  EPASF  (Reference  20). 

vii.  Two-Dimensional  Radial  Transport  Model  LTIRD 
(Reference  16). 

b.  Semi-analytical  Model 

I.  Two-Dimensional  Complex  Velocity  Potential  Model  RESSQ 
(Reference  16). 

c.  Numerical  Models 

i.  Method  of  Characteristics  Model  MOC  (References  21,  22 
and  23). 

ii.  Random  Walk  Solute  Transport  Model  RWALK 
(Reference  24). 

Several  investigators  have  compared  the  performance  of  numerical  codes  to 
analytical  solutions,  benchmark  data  sets  and  real  site  applications  (References  25, 
26,  27,  28  and  29).  The  algorithm  for  choosing  among  the  numerical  codes  is  based 
in  part  on  such  comparisons.  Another  version  of  the  algorithm  is  under  development, 
capable  of  selection  of  transport  models  used  to  predict  the  effectiveness  of 
alternative  remediation  schemes,  optimizing  for  cost/effectiveness.  In  essence,  the 
first  algorithm  suggests  a  model  or  models  for  the  initial  transport  prediction:  the 
second  will  provide  guidance  on  the  remediation  method,  and  this  may  in  turn  require 
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selection  of  another  transport  code. 

The  algorithm  is  outlined  in  Figure  5.  The  user  responds  to  screen  queries 
about  whether  analytical  solutions  are  known  to  be  appropriate  or  inappropriate  (in  the 
latter  case,  whether  the  region  modeled  is  homogeneous  or  heterogeneous).  If  the 
complexities  require  a  numerical  model,  the  algorithm  then  jumps  to  that  branch  to 
select  between  the  two  available  numerical  codes.  The  flow  charts  identify  the  model 
recommended  as  a  result  of  certain  user  responses:  whether  the  subsurface  waste 
disposal  method  is  a  landfill,  a  wastewater  lagoon  or  spray  irrigation;  whether  the  flow 
is  radial  or  not,  whether  the  Dupuit  approximation  is  valid  or  not;  whether  single  or 
multiple  sources  are  involved;  whether  full  penetration  analysis  is  adequate;  whether 
regional  flow  is  important  or  not.  For  example,  the  algorithm  checks  if  a  particular 
solution  applies:  if  the  user  responds  in  the  affirmative  that  flow  in  the  region  is 
strongly  affected  by  pumping  wells,  then  semi-analytical  (complex  velocity  potential) 
methods  or  complete  numerical  methods  would  be  indicated  as  more  appropriate  than 
models  based  upon  analytical  solutions.  In  the  case  of  selecting  a  numerical  model, 
the  user  is  prompted  to  respond  to  queries  about  grid  size,  longitudinal  and  transverse 
dispersion,  whether  the  flow  is  parallel  to  the  grid  axes,  whether  storativity  is 
significant,  and  whether  any  part  of  the  aquifer  changes  from  confined  to  unconfined 
flow  or  vice  versa. 
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Vigura  5.  CHOICE  Algorithm  Flow  Chart 
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rigur«  5.  choice  Algorithm  Flow  Chart  (cont'd) 
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rigurtt  S.  CHOICE  Algorithm  Flow  Chart  (cont'd) 


21 


22 


rigur*  S.  CHOICE  Algorltha  Flow  Chart  (cont'd) 
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Figur*  S.  CHOICE  Algoritha  Flow  Chart  (cont'd) 
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2. 


LeGrand  Method 


For  the  purposes  of  preliminary  site  analysis,  we  have  included  in  the 
System  the  standardized  evaluation  methodology  developed  by  LeGrand 
(Reference  30).  The  LeGrand  method  is  not  a  contaminant  transport  model,  but  is  a 
tool  for  the  preliminary  assessment  of  a  site  given  certain  commonly  available 
information  on  the  site  and  its  environment.  The  methodology  offers  a  concise 
screening  mechanism  for  evaluating  the  contamination  potential  of  waste  sites,  as  well 
as  a  management  control  procedure  useful  during  planning  and  operational  stages  of 
contaminant  handling.  It  is  not  as  comprehensive  as  the  Air  Force  Automated  Defense 
Priority  Model  (ADPM). 

The  LeGrand  method  focuses  on  weighting  appropriately  the  key  characteristics 
of  a  site,  in  a  standardized  manner  to  form  a  preliminary  evaluation  of  contamination 
potential.  Each  key  characteristic  is  assigned  a  numerical  value.  The  method  relies 
on  the  quantification  of  certain  parameters,  evaluated  in  a  logical  sequence,  with  the 
results  presented  in  a  standardized  form.  The  numerical  rating  system  is  divided  into 
ten  steps  within  four  stages: 

e  Standard  hydrogeological  description  of  the  site.  (Stage  I) 

e  Determination  of  how  serious  the  hazard  potential  is  by  identifying  the 
degree  of  aquifer  sensitivity  and  the  degree  of  contaminant  severity. 
(Stage  II) 

•  Description  of  the  relative  probability  of  contamination  by  comparing  the 
site's  numerical  value  with  a  standard  value  that  is  derived  from 
consideration  of  both  aquifer  sensitivity  and  contaminant  severity. 
(Stage  III) 

e  Reassessment  of  the  site,  with  consideration  given  to  engineering 
modifications.  (Stage  IV) 


Stage  I:  Numerical  Description  of  Site  HvdroaeoloQv 

The  numerical  rating  system  first  provides  a  description  of  the  site 
hydrogeology.  This  is  based  primarily  on  four  key  characteristics: 

a.  Distance  on  the  ground  from  a  source  of  contamination  to  the  nearest 
well,  surface  stream  or  property  boundary; 

b.  Depth  of  the  water  table  below  the  waste  or  contamination  source; 
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c.  Approximate  slope  of  the  water  table:  and 

d.  Character  of  the  soil  materials  through  which  the  contaminant  is  likely  to 
pass,  expressed  in  terms  of  permeability  and  sorption.  The  method  to  describe  site 
hydrogeology  assigns  a  ”0”  value  to  the  most  favorable  setting  of  each 
hydrogeological  factor,  and  a  ”9”  value  to  the  least  favorable.  Intermediate  numerical 
values  are  defined  by  interpolating  conditions  between  the  two  extremes.  For  each 
site  the  estimated  numerical  or  point  value  for  each  of  the  four  factors  is  added  and 
the  total  expressed  as  a  number  between  0  and  32  that  characterizes  the  site. 

This  stage  consists  of  seven  steps.  The  first  four  steps  involve  the  recording 
of  estimated  values  for  each  of  the  four  key  hydrogeological  parameters.  These  four 
steps  represent  the  core  of  the  description  method* 

1 .  Distance.  A  great  distance  between  a  contamination  source  and 
a  water  supply  (or  perimeter  of  regulatory  compliance)  is  generally  considered  a 
favorable  factor,  especially  in  loose,  granular  materials  with  some  sorption  capacity. 
It  is  a  less  significant  factor  in  fractured  materials  or  cavernous  rock,  where 
contamination  is  likely  to  reach  gr'^ater  distances.  To  distinguish  one  situation  from 
the  other  the  identifying  letter  T  is  introduced  in  Step  6. 

2.  Depth  to  the  water  table.  Many  contaminants  attenuate  in  the 
unsaturated  zone  above  the  water  table.  Therefore,  it  is  generally  beneficial  to  have 
a  deep  water  table.  The  depth  to  the  water  table  is  defined  here  as  the  distance  from 
the  ground  surface  to  the  surface  of  water  standing  in  an  unpumped  well  (pumped 
wells  will  lower  the  water  table  locally).  The  points  assigned  increase  with  shorter 
distance  to  the  water  table.  The  scale  is  not  a  simple  arithmetic  progression  in  order 
to  account  for  the  increased  sensitivity  at  smaller  depths  to  the  water  table. 

3.  Water  table  gradient.  This  factor  establishes  whether  the 
contaminant  is  moving  toward  or  away  from  the  water  supply  (or  perimeter  of 
interest).  Since  the  exact  gradient  is  difficult  to  estimate  unless  current  water  table 
maps  are  available,  only  5  points  are  spread  across  this  scale. 

4.  Permeability  and  sorption.  The  scale  for  the  permeability-sorption 
factor  extends  from  0,  representing  the  low  permeability  and  high  sorption 
characteristic  of  clay,  to  9,  representing  the  high  permeability  and  low  sorption  of 
clean  gravel.  The  sites  are  also  classified  using  a  letter  qualifier.  The  letter  attached 
to  each  digit  in  the  matrix  of  permeabiiity-sorption  is  for  specific  identification  of  the 
characteristics  of  a  site. 

5.  Steps  5  and  6  provide  for  the  addition  of  letter  identifiers  that 
identify  special  features  with  respect  to  the  site. 
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6.  In  Step  7  the  separate  ratings  and  identifying  letter  suffixes  are 
recorded  and  the  values  are  added  to  achieve  a  total  description. 

A  site  description  from  the  LeGrand  method  is  a  compound  of  four  separate 
digits,  representing  values  from  the  first  four  steps,  and  four  or  more  letters  derived 
from  Steps  4,  5  and  6.  The  first  letter  identifier  is  derived  from  the  permeability 
sorption  matrix  (Step  4).  The  second  letter  (A,  B  or  C)  assigns  a  level  of  confidence 
to  the  overall  values  derived  from  Steps  1  to  4  (Step  5).  The  third  letter  indicates 
whether  the  distance  from  a  contamination  source  is  to  a  well,  spring,  perennial 
stream  or  specified  boundary.  The  fourth  letter  identifier  is  selected  from  the  most 
appropriate  characteristic  listed  with  the  letter  identifiers  in  Step  6;  an  additional  letter 
identifier  from  this  list  may  also  be  added. 

Step  7  completes  the  site  numerical  description.  It  is  accomplished  by  adding 
the  separate  point  values  determined  in  the  first  four  steps  and  writing  the  sum  with 
the  appropriate  number  value  and  letter  suffixes.  At  this  stage  the  site  can  be  rated 
in  terms  of  relative  hydrogeological  conditions,  but  not  necessarily  with  respect  to  the 
possibility  of  contamination.  The  site  is  assigned  a  grade  based  on  its  critical 
hydrogeological  parameters  as  assessed  in  Steps  1-6. 


Stage  II;  Evaluation  of  Degree  of  Seriousness 

The  total  hazard  potential  has  two  components:  degree  of  seriousness  and 
probability  of  contamination.  The  degree  of  seriousness  is  somewhat  independent  of 
the  site  description  obtained  in  Stage  I,  but  is  an  essential  part  of  any  groundwater 
contamination  evaluation.  The  analysis  performed  at  this  level  is  made  with  a  matrix 
considering  the  sensitivity  of  the  aquifer  to  contamination  and  the  severity  of  the 
contaminant.  These  factors  are  defined  as  follows: 

a.  Aouifer  sensitivity.  This  term  is  used  to  indicate  likelihood  and  degree 
of  groundwater  contamination  at  a  given  site.  The  aquifer's  areal  extent  and 
importance  as  a  ground  water  source  are  also  considered  in  this  definition.  For  use 
in  this  rating  system,  aquifers  are  divided  into  three  categories:  sensitive,  moderately 
sensitive  and  insensitive.  Permeability  is  the  key  factor  in  considering  aquifer 
sensitivity.  Additional  factors  are  the  thickness  of  the  saturated  part  of  the  aquifer 
and  the  quality  of  water  with  respect  to  its  acceptability  for  use. 

b.  Contaminant  severity.  This  term  includes  qualitative  weighting  of 
toxicity,  concentration  and  volume,  mobility  in  the  groundwater  and  persistence.  The 
contaminant  severity  scale  ranges  from  the  effluent  of  a  single  septic  tank  at  the  low 
end  to  large  volumes  of  high-level  radioactive  wastes  at  the  high  end. 

The  overall  degree  of  seriousness  is  determined  by  the  intersection  between 
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"contaminant  severity"  and  "aquifer  sensitivity."  The  overall  degree  of  seriousness 
is  divided  into  nine  categories,  ranging  from  "relatively  small”  to  "extremely  high." 


Stage  III:  Evaluation  of  Probability  of  Contamination 

The  matrix  of  contaminant  severity  and  aquifer  sensitivity  is  used  again  in 
Stage  III  to  grade  a  situation  more  specifically.  For  each  intersection  of  these  two 
parameters  a  standard  situation  rating  is  defined.  The  standard  situation  rating 
represents  the  approximate  numerical  value  which  a  site's  description  should  not 
exceed  to  prevent  serious  contamination  of  groundwater.  These  PAR  values  (Pro¬ 
tection  of  Aquifer  Ratings)  were  derived  from  extensive  studies  of  a  large  number  of 
situations  of  varying  aquifer  sensitivities  and  contaminant  severities,  including  cases 
where  contamination  has  or  has  not  occurred.  These  values,  empirical  in  nature,  are 
then  compared  with  the  hydrogeological  numerical  grade  for  the  site  obtained  in 
Stage  I.  Based  on  this  comparison  we  obtain  a  situation  rating,  from  which  the 
probability  of  contamination,  degree  of  acceptability  and  situation  grade  are  derived 
for  the  site  under  consideration.  Stage  III  consists  of  Steps  8  and  9.  Step  8  evaluates 
the  probability  of  contamination  and  degree  of  acceptability  for  a  natural  site,  and 
Step  9  performs  this  evaluation  for  a  modified  site. 


Stage  IV;  EnoineerinQ  Improvements  and  Final  Acceptance 

Generally,  the  areas  around  potential  contaminant  sites  are  heavily  influenced 
by  human  modifications  which  result  in  changes  in  the  subsurface  hydrologic  regime. 
Engineering  modifications  to  limit  contamination  are  common  and  should  be 
evaluated.  Stage  IV  provides  a  means  to  rate  sites  that  are  modified  by  human  action. 
These  changes  result  in  aquifer  sensitivity  and/or  contaminant  severity  modifications. 
In  Stage  IV  the  effect  of  the  modified  properties  is  reassessed  and  the  new  PAR  values 
are  evaluated.  This  results  in  a  new  situation  rating,  giving  new  probability  of 
contamination,  degree  of  acceptability  and  situation  grade.  This  stage  provides  a 
simple  method  to  predict  the  impact  of  human  action  on  the  contamination  potential 
of  a  site,  and  establish  a  preliminary  evaluation  of  the  cost/benefit  ratio. 


3.  ODAST 

The  model,  ODAST,  documented  in  Javandel  et  al.  (Reference  16), 
provides  an  analytical  solution  of  the  one-dimensional  convective-dispersive  transport 
equation,  and  is  adapted  from  the  set  of  one-dimensional  solutions  published  by  Van 
Genuchten  and  Alves  (Reference  31). 

For  the  derivation  of  the  model,  consider  a  one-dimensional  system 


29 


approximated  as  an  infinitely  long  homogeneous  isotropic  porous  medium  with  a 
steady  state  uniform  flow,  with  seepage  velocity  v.  A  dissolved  constituent  is  injected 
at  one  end  of  the  model  for  a  period  of  time  to,  such  that  the  input  varies  as  an 
exponential  function  of  time.  The  governing  equation  for  this  situation,  in  a  confined 
homogeneous  isotropic  aquifer  (ignoring  the  dependence  of  dispersion  on  space  and 
time),  is: 
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dx  dt 


where  D  is  the  dispersion  coefficient,  lambda  (A)  is  a  radioactive  decay  constant,  and 
R  is  the  retr'd‘:tion  coefficient.  This  equation  is  solved  subject  to  the  following  initial 
and  boundary  conditions: 

C{x,t)=0  /  =  0 

=0  jr  =  00 


dx 


vC 


U.0  =  ^  nt) 


t)  =  Cq  exp  (  -at)  0  <  /  <  /o 

/(/)  =  0  f  >  /o 

where  Cq  and  a  are  constants. 

Using  the  Laplace  transform  technique.  Van  Genuchten  and  Alves 
(Reference  31)  solved  this  set  of  equations,  obtaining: 

d7(jr, /)  =  /l(x,t)  0  ^t  ^to 


c(x,t)  =A(x,t)  -Ai(x,t-/o)  exp(-afo)  t  > 


where: 


A{x,t)  =  CJ  exp( -a)  v4,( jf. /)  a^X 


A{x,  t)  =  Q  exp(  -at)  A^{x,t)  5  =  X 


t)  = 


V  *  U 
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The  model  can  then  be  adapted  for  Monte  Carlo  simulation  by  treating  D  and 
V  as  random  variables.  The  velocity  {v)  is  not,  however,  directly  generated:  it  is 
assumed  to  result,  via  Darcy's  law,  from  hydraulic  gradient,  porosity  and  hydraulic 
conductivity.  Hydraulic  conductivity  may  be  generated  directly  for  the  simulation  via 
a  log-normal  probability  distribution.  However,  this  method  does  not  take  into  account 
the  known  relationship  between  the  hydraulic  conductivity  and  other  parameters. 
Therefore,  the  option  is  also  provided  in  the  code  to  generate  a  distribution  of 
hydraulic  conductivities  from  the  underlying  variables  of  mean  particle  size  and 
porosity,  from  the  Karmen-Cozeny  relationship.  This  is  the  approach  most  appropriate 
for  preliminary  analysis  of  contamination  risk  when  there  is  considerable  uncertainty 
regarding  the  hydrogeology.  The  underlying  variables  of  gradient  and  mean  particle 
size  are  generated  for  Monte  Carlo  simulation  in  the  manner  described  for  the  model 
EPAGW. 
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4.  TDAST 


The  model  TDAST  provides  an  analytical  solution  to  steady  state  flow  In 
a  two-dimensional  Cartesian  coordinate  system,  and  is  documented  in  Javandel  et  al. 
(Reference  16)  If  the  flow  is  coincident  with  the  x  axis,  and  the  longitudinal  and 
transverse  components  of  the  dispersion  tensor  are  assumed  independent  of  position 
and  designated  by  D^^  and  Dj,  the  general  governing  equation  for  a  confined, 
homogeneous,  isotropic  aquifer  can  be  written  as: 


If 

dx  dt 


where  lambda  [A]  represents  decay  and  R  is  the  retardation  coefficient. 


For  a  particular  solution,  we  first  assume  that  the  medium  is  initially  free  of  the 
solute  and  that,  at  a  certain  time,  a  strip  type  source  of  length  2a,  orthogonal  to  the 
flow  direction,  is  introduced  along  the  y  axis.  If  the  source  concentration  diminishes 
exponentially  with  time,  the  initial  and  boundary  conditions  are: 


C{o,y,  t)  =  Cq 
o,  y,t)  =0 


lim  dC 

y-*±co 


=  0 


lim  dC 


=  0 


-a  <y  <  a 
\y\  >  a 


Where  the  source  "strip"  is  arranged  orthogonal  to  the  direction  of  flow,  an 
analytical  solution  is  presented  by  Cleary  and  Ungs  (Reference  32)  as: 


This  model  can  be  adapted  for  Monte  Carlo  simulation  in  a  manner  analogous 
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to  ODAST,  except  thet  here  the  variable  Dj  must  also  be  generated. 
S.  PLUM2D 


The  model  we  refer  to  as  PLUM20.  originally  PLUME  2D,  was  developed 
to  provide  an  analytical  solution  for  multiple,  fully  penetrating  point  sources  in  a 
homogeneous,  nonleaky  confined  aquifer  with  uniform  regional  flow,  with  the  flow 
coincident  with  the  x  axis.  The  solution  is  documented  in  van  der  Heijde 
(Reference  17).  Sources  are  conceived  as  fully  penetrating  wells.  In  this  case,  the 
general  governing  equation  for  a  solute  subject  to  radioactive  decay  and  adsorption 
described  by  a  linear,  equilibrium  relationship  can  be  written  as: 


R  -  D  +  D 


-  p 

dx 


Xf%c 


where  the  dispersion  coefficients  may  be  written  as: 

Dx  =  a V  \  Dy^OyV 

where  — 
n 

where  v  is  the  Darcy  velocity,  n  is  the  porosity,  o,  and  Oy  are  the  dispersivity  In  the  x 
and  K  direction,  is  the  retardation  coefficient,  and  A  =  In2/T,  with  r  being  the  half- 
life  time. 


For  an  infinite  two-dimensional  porous  medium  in  the  x,y  plane,  with  a  mass  per 
unit  length  in  the  ^-direction,  Mj,  instantaneously  injected  along  the  vertical  z-axis,  the 
solution  is  given  by: 


C{x,y,t) 


4nnt  y]  UjfUy 


exp 


The  solution  may  be  written  by  analogy  to  Hantush's  leaky-well  function, 
Wfu,r/BJ,  following  Wilson  and  Miller  (Reference  33),  as: 


C{x,y.t) 


AnnfqUy 
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where  verieble  B  forms  the  mixing  scale  equivalent  of  the  Hantush  leakage  factor, 
which  accounts  for  the  effects  of  dispersion  in  the  solute  transport  model. 


Because  of  the  nature  of  the  solution,  various  sources  may  now  be 
superimposed  to  estimate  the  result  of  contamination  from  multiple  point  sources, 
which  may  have  been  operational  for  varying  lengths  of  time.  The  solution  used  in  the 
program  is  obtained  by  use  of  an  accurate  analytical  approximation  of  the  Hantush 
ieaky-weil  function,  given  by  Walton  (Reference  34).  The  model  may  then  be  adapted 
for  Monte  Carlo  simulation  in  a  manner  similar  to  that  for  other  analytical  solutions 
derived  here. 

6.  DUPVG 

The  model  denoted  as  DUPVG  in  the  system  was  originally  developed  to 
address  the  problem  of  solute  transport  in  an  unconfined  aquifer,  where  the  water 
table  surface  may  be  moving  in  response  to  input  from  a  source  (References  18  and 
35).  The  equations  governing  flow  and  transport  in  the  saturated  zone  of  unconfined 
aquifers  are  subject  to  the  non-linear  moving  free-surface  boundary.  Where  the 
movement  of  this  free  surface  is  a  significant  component  of  the  flow  regime  analytical 
solutions  based  on  the  assumption  of  a  confined  aquifer  are  no  longer  appropriate. 
The  model  of  Volker  and  Guvanasen  provides  an  approximate  analytical  solution  for 
such  cases. 

Consider  a  case  in  which  an  infinitely  long  (in  y  direction)  unconfined  aquifer  is 
being  recharged  through  an  infinitely  long  recharge  basin  of  width  2L.  A  symmetrical 
setting  is  assumed,  and  the  aquifer  is  bounded  by  a  drain  of  constant  head  on  each 
side  (x  direction)  at  distance  B  from  the  center  of  the  recharge  basin.  The  basin  is 
recharging  the  aquifer  at  rate  P^.  The  medium  is  assumed  homogeneous  and  uniform 


with  effective  porosity  6.,  hydreulic  conductivity  K  and  initial  saturated  depth  Fiow 
is  assumed  to  arise  only  as  a  result  of  the  recharge  (/?),  which  results  in  a  spatially 
varying  head-rise  of  s. 


By  symmetry,  the  solution  need  be  sought  only  in  the  x  >  0  segment.  The 
procedure  adopted  by  the  authors  is  first  to  solve  the  flow  equations,  governing  the 
velocity  distribution,  then  appiy  these  to  the  transport  equation.  The  rise  of  the  free 


surface  above  the  initial  saturated  depth, . 
Forcheimer  assumption.  The  resulting 
assumption  that  s/a^  <<  1,  which  will 
equation  governing  flow  is  then: 

dx^ 

subject  to  boundary  conditions: 

5=0. 
ds!  dx  =  0, 

5=0. 

/?=  Po, 

/?=  0. 


t,  is  first  obtained  by  applying  the  Oupuit- 
equation  is  then  linearized  under  the 
be  true  for  small  infiltration  rates.  The 

_  ^ 
dt 


/  =  0.  Orsx 

X  =  0 

X  ^  B 

0  ^  X  s  ^  Z. 

L  ^x  ^B 


A  solution  to  this  equation  can  be  obtained  through  the  eigen-function 
expansion  method,  defining  the  position  of  the  upper  boundary  of  the  saturated  flow 
domain.  However,  it  should  be  noted  that  the  Dupuit-Forcheimer  assumption  implies 
a  horizontal  flow  beneath  the  basin,  where  a  vertical  flow  component  may  actually 
predominate.  To  obtain  a  more  accurate  representation  of  the  flow  field  the  authors 
first  assume  that  the  flow  pattern  at  any  time  t  can  be  described  as  a  function  of  the 
steady  state  velocity: 

Ui  (X.  2,  t)  =  U,  {X,  2.0^  f(t) 


where  f(t)  is  the  scale  function  dependent  on  time.  The  velocity  distribution  is  then 
sought  by  transforming  the  free-surface  flow  domain  to  rectangular  coordinates: 


X  =  x. 


Zs 


5  +  5o 


Given  the  assumption  that  5/ao  <  <  /,  it  can  be  further  assumed  that  the 
unsteady  free  surface  can  be  approximately  described  by  a  streamline,  that  streamline 
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and  equipotential  functions  change  tittle  with  time,  and  that  further  away  from  the 
source  the  velocity  is  essentially  horizontal  and  its  spatial  variation  is  negligible.  Based 
on  these  assumptions,  only  the  steady-state  velocity  pattern  in  the  transformed 
domain  need  be  sought,  and  the  downstream  end  can  be  extended  to  infinity.  The 
flow  equation  and  associated  boundary  conditions  are  simplified  to: 


Solutions  are  then  obtained  by  the  Schwarz-Christoffel  transformation  method 
in  terms  of  the  relationships  between  X,  Z,  and  the  equipotential  and  stream  funcliwi., 
yielding: 

jf.  ^  In  (v^  .  ./rTTET) 


z  - 


—  COS'^ 
Tt 


Y, 


cosh  {^) 


0  =  .P  In  {  /v  +  >/(  1  +v)  } 


where 

Yi  =  cosh  (  ^)  cos  ( ■^)  +  €2 

Y2  =  €,  sinh  (  sin  ( 

5  =  -^  {  (  yI  +  Y1-I)  +  /(yI  Y?-1)^  4y2  +  4Yi  } 
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c, 

X,  =  ±  cosh  {^)  cos  {^)  -  Cg/C, 

*  1  ®0  *0 

=  ^  sinh  sh  {^) 

*  1  «0  ®0 

v  =  -g\lX|  +  Ai-1)  +  a|  +  Xi-2+  4>|  } 


The  steady  state  velocity  along  streamlines  is  then  given  by: 
V  (0.  /-►oJ  * 


O 

€l«0 


( Yi  ^  yI-I)^  *  4y2 


2  11/4 


^  cosh  2(  S//7  2(  4  sAlh  2( 


By  assumption: 


u^{f)  = 

f(t)  =1  -  exp  I  ■ 


and 

-Ky^ 


The  transformation  of  the  domain  is  then  applied  to  the  transport  equation,  with 
the  additional  transform  of  time: 

T  =  no  <Ao 


With  the  assumptions  that  s/a^  <<  1.  the  slope  of  the  free  surface  is  small,  the  rate 
of  rise  of  the  free  surface  is  very  small,  and  there  is  no  dispersion  across  streamlines, 
the  transport  equation  can  be  expressed  in  the  curvilinear  equipotential-streamline 
coordinate  system  as: 


dC 

dr 


+ 


2  ^ 

dP 
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For  use  in  the  Advisory  System  we  follow  Solution  SI  as  given  by  Guvanasen  and 
Vollcer  (Reference  35).  The  reduced  governing  transport  equation  is  solved  subject  to 
the  following  conditions: 

(1)  The  boundary  condition  upstream  is  equivalent  to: 

c{0  =  o.T)  =  a 

which  implies  a  concentration  gradient  across  the  boundary  approaching  zero. 

(2)  At  the  vertical  downstream  end,  with  0  <«  0^: 

dC  _  C 

^  ‘  TDl 


which  is  an  approximation  for  lower  values  of  C. 


(3)  The  velocity  along  streamlines  at  steady  state  is  assumed  to  be 
uniform  and  equivalent  to: 


which  assumes  that  can  be  taken  outside  the  differentiation,  and  is  thus 
independent  of  position. 

(4)  Initial  condition: 

C{0.T  -0)  =0. 


By  Laplace  transform  methods,  Volker  and  Guvanasen  (Reference  1 8)  show  that 
the  solution  then  becomes: 


(  -0  +  2/70^ 

2D, 


.  (  0-2/700 -ir*T) 

erfc  AZ - ZE. — IL_L  4  exp 
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0-2^0g  \  .  (<p-2n<Pg* 


2Zr  j 


Noting  that  B,  the  distance  to  the  constant  head  boundary,  goes  to  infinity,  a  good 
approximation  to  the  equation  above  is: 
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Brfc 


erfc 


2i^»v^  J 


This  approximation  will  not  be  particularly  good  when  B  is  small  relative  to  the 
distance  from  the  source.  However,  this  approximation  is  only  involved  in  the  solution 
of  the  transport  equation,  and  not  in  the  solution  of  the  velocity  equation.  Thus  the 
error  should  not  affect  the  computed  average  position  of  the  front. 

7.  EPAGW,  Modified  EPA  Model  for  Monte  Carlo  Analysis  of  impact  On 

Groundwater 

This  groundwater  model  accounts  for  most  of  the  major  physical  and 
chemical  processes  that  affect  movement  and  transtormations  of  chemicals  in  simple, 
homogeneous  and  isotropic  porous  media  under  steady  flow  conditions 
(Reference  20).  The  mechanisms  considered  include  advection,  hydrodynamic 
dispersion  in  the  longitudinal,  lateral  and  vertical  dimensions,  adsorption  and  chemical 
degradation.  Major  assumptions  made  by  the  transport  model  are: 

a.  groundwater  velocity  is  uniform,  one-dimensional  and  steady-state, 
in  a  saturated  aquifer; 

b.  the  porous  medium  is  homogeneous  and  isotropic; 

c.  mass  transport  is  also  in  a  steady  state; 

d.  the  aquifer  is  semi-infinitely  large  in  the  direction  of  groundwater 
flow  and  infinitely  large  In  the  transverse  direction; 

e.  the  contaminant  source  is  sufficiently  large  in  mass  such  that  the 
down-gradient  concentration  will  be  maintained  once  it  is  reached; 

f.  the  distribution  of  contaminant  concentration  at  the  source 
boundary  is  Gaussian; 

g.  degradation  of  chemicals  is  caused  only  by  hydrolysis; 

h.  equilibrium,  reversible  speciation,  and  sorption  are  appropriate  -- 
in  order  to  utilize  an  equilibrium  partition  coefficient. 

The  advection-dispersion  equation  for  the  transport  of  a  nonconservative 
contaminant  in  an  adsorbing  homogeneous  and  isotropic  porous  medium  with  fully- 
saturated  flow  may  be  written  as  follows: 


39 


where 


3jr* 

V 


dx 


dC 


D„±£  - 
^  dz^ 


^  ^  ^  ^  xc*  tc 


x,y,z 

Spatial  coordinates  in  the  longitudinal,  lateral  and  vertical 
directions,  respectively,  fUi 

C 

= 

dissolved  concentration  of  chemical,  [MA\ 

retarded  dispersion  coefficients  in  the  x,y  and  z  directions, 
respectively,  (L*/T)', 

V 

groundwater  seepage  velocity,  assumed  to  be  in  the  x 
direction,  (Zy7); 

= 

retardation  factor; 

t 

as 

elapsed  time,  (7); 

e 

as 

volumetric  water  content  of  the  porous  medium; 

1 

= 

net  recharge  due  to  precipitation,  (7*'). 

The  retardation  factor  and  the  effective  reaction  rate  constant  are  defined  as 

follows: 

-  1  •  ^ 

^  0  +  X2  p^ 

0  *  P,  K 

where  p^, 

s 

bulk  density  of  the  porous  medium, 

= 

distribution  coefficient, 

0 

as 

volumetric  water  content; 

s 

rate  coefficient  for  dissolved  phase,  { f/7); 

^2 

=s 

rate  coefficient  for  sorbed  phase,  (7/7). 

The  three-dimensional  region  of  interest  is  regarded  as  semi-infinite  in  the  x- 
direction  (0  ^  x  <  <»),  Infinite  in  the  y-directlon  (-»  <  y  <  <»),  and  finite  In  the  z- 
direction  (0  ^  z  ^  B),  where  B  is  equal  to  the  saturated  zone  thickness.  The  above 
partial  differential  equation  is  solved  analytically.  The  solution  treats  the  source 
concentration  as  a  Gaussian  distribution  in  the  lateral  direction  (along  the  y-axis 
corresponding  to  the  leading,  down-gradient  edge  of  the  unit),  and  a  uniform 
distribution  over  the  vertical  mixing  or  leachate  penetration  depth,  H.  The  maximum 
dissolved  concentration  of  contaminant,  C„,  occurs  at  the  center  of  the  Gaussian 
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distribution,  y^.  This  distribution  is  defined  by  its  standard  deviation,  a,  which  is 
measured  in  terms  of  distance  (in  meters)  and  is  related  to  the  width  of  the  disposal 
unit.  The  initial  and  boundary  conditions  necessary  to  solve  the  equation  may  be 
written  as  follows: 

C(x,y,z,0)  =0 

C{o,y.z,t)  =  Q  exp ( -/*/ 2o2)  U{  Zi 
C{x,  «nz,t)  =0 
C{x,  -on  z,t)  =0 
C{c>{y,z.  /)  =  0 

^Uy.o.n-0 

^(x.y.st)  =0 

where  Co  is  the  peak  concentration  at  the  source,  a  is  the  standard  deviation  of  the 
Gaussian  distribution  centered  at  x  =  y  =  0.  and  U(z)  is  a  unit  step  function 

U(z)  =1  .  Hi^z^fi 

U{z)  =0  ,  z>  or  z> 


A  steady-state  analytical  solution  can  be  obtained  by  direct  solution  of  the 
steady-state  version  of  the  governing  differential  equation,  which  implies  removal  of 
the  time-derivative  term.  Details  of  the  solution  procedure  are  contained  in  a  U.S.  EPA 
report  (Reference  36).  The  solution  is: 

Cp(x,y,z)  =^Cf[x,/i  ^bCp{x,y,z) 


where 


/  f 


Cf{x,y)  =5'  J 


exp  (-y|/2ag) 
yj  DJ  Dy 


m!  +  ^  ^  Yp-y) 


4 


41 


arKi 


in  which  is  a  dummy  variabie  of  integration^  and  f  is  a  constant: 


and  K^{’)  is  the  modified  Bessei  function  of  the  second  kind  and  of  first  order.  The 
solution  given  assumes  that  the  source  extends  from  the  top  of  the  aquifer  through 
the  thickness  H,  so  that  W,  =0  and 

To  obtain  the  steady-state  concentration  distribution  along  the  x  (flow)  axis,  the 
solution  can  be  reduced  to: 

Cp  (x.0,0)  =  -g  Cfi  X,  0)  +  ACp  ( X,  0.  0) 

where: 

c,(x.O) 
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El-Kadi  et  al.  (Reference  37)  applied  this  model  to  six  waste  disposai  sites  in 
different  regions  of  the  United  States.  The  modei  was  tested  against  a  three- 
dimensional  finite  element  model  (CFEST)  and  was  found  to  generally  underestimate 
the  contaminant  concentrations:  they  recommend  using  a  factor  of  safety  of  one  order 
of  magnitude  for  a  conservative  analysis.  It  should  be  emphasized  that  the  "exact” 
results  were  taken  to  be  those  produced  by  the  numerical  model.  Overestimation  of 
concentrations  resulted  when  the  modei  was  applied  to  a  site  where  radial  flow 
existed.  With  appropriate  choice  of  parameters,  analytical  solutions  are  especially 
suitable  for  management  purposes  when  combined  with  sensitivity  analysis  and/or 
uncertainty  analysis  (e.g.,  Monte  Carlo  simulation). 

Input  data  for  all  model  parameters  must  be  identified  to  use  the  EPAGW 
model.  In  general,  however,  the  behavior  of  a  specific  constituent  in  the  environment 
is  highly  dependent  on  both  the  environmental  setting  and  the  properties  of  the 
constituent.  The  assignment  of  specific  values  to  describe  the  behavior  of  the 
modeled  system  is  further  complicated  by  uncertainty  about  how  to  specify  a  single 
value  for  each  model  parameter  which  represents  a  "worst  case”  (obtained  from  the 
steady-state  solution  along  the  x-axis). 

As  an  alternative  to  identifying  reasonable  worst-case  values  for  each  model, 
a  Monte  Carlo  simulation  technique  is  used.  It  involves  a  large  number  of  computer 
runs  with  values  for  each  input  parameter  drawn  from  data  sets  describing  ranges  of 
possible  values  and  the  distribution  of  values  within  the  range.  Where  parameters  are 
correlated,  and  therefore  dependent,  the  relationships  are  properly  defined  in  the 
Monte  Carlo  routine.  The  Monte  Carlo  process  in  EPAGW  (and  in  general)  proceeds 
as  follows: 

(1)  Values  from  each  input  distribution  are  selected  at  random. 

(2)  The  model  is  run  with  a  set  of  randomly-generated  parameters  to  give  a 
set  of  output  variables. 

(3)  The  input  selection  and  computation  steps  are  repeated  a  large  number 
of  times  (1000  to  5000)  to  produce  a  well-defined  output  distribution. 

(4)  The  output  values  are  analyzed  for  presentation  as  distribution. 
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The  groundwater  model  parameters  and  input  data  requirements  include  the 
following:  groundwater  velocity,  porosity  of  the  saturated  media,  dispersivity  of  the 
aquifer,  distance  to  the  measurement  point,  standard  deviation  of  the  Gaussian 
source,  penetration  depth  of  leachate  into  the  aquifer,  thickness  of  aquifer,  fraction 
of  organic  carbon  content  of  the  soil,  pH  and  temperature  of  groundwater,  acid, 
neutral  and  base  hydr^'lysis  rates.  Relationships  between  these  environmental 
parameters  must  be  determined  in  order  to  apply  the  Monte  Carlo  analysis  properly. 
All  the  parameters  and  input  data  previously  mentioned  are  in  some  extent  linked  or 
dependent  on  at  least  one  of  the  others.  In  some  cases  an  independent  seed 
distribution  can  be  generated  to  which  other  variables  are  correlated.  Such  is  the  case 
of  temperature,  which  is  not  a  function  of  any  of  the  other  parameters,  but  influences 
the  values  of  some  of  them  (e.g.,  hydrolysis  rate  coefficients).  Dependent  data  sets 
can  be  developed  as  empirical,  joint  or  multivariate  distributions,  theoretical 
distributions,  or  from  functional  deoendenc'es  among  the  variables  and  parameters. 
The  parameters  and  variables  to  be  generated  independently  are  as  follows: 

a.  Thickness  of  the  saturated  zone,  B 


b.  Fractional  organic  carbon  content,  FOC,  The  fractional  organic 
carbon  content  is  used  to  determine  the  distribution  coefficient,  Kq,  from  the  following 
relationship: 


/CDo  =  (foe)  (Koc) 


where 


distribution  coefficient  normalized  to  organic  carbon. 


c.  Groundwater  pH,  assumed  to  be  independent  of  contaminant 
concentration  and  temperature. 


d.  Groundwater  temperature,  T. 

e.  Leachate  penetration  into  the  saturated  zone,  H,  probably  related 
to  the  relative  differences  in  the  leachate  and  groundwater  velocities.  A  simple. 
Independent,  uniform  distribution  ranging  from  a  fixed  minimum  to  a  fixed  maximum 
is  used  to  represent  this  variable. 


f.  Net  recharge,  /.  It  is  the  amount  of  water  that  enters  an  aquifer 
system.  Since  the  groundwater  model  assumes  that  the  porous  media  is  uniform,  the 
effect  of  recharge  causes  the  groundwater  to  rise  and  fall  uniformly.  There  is  thus  no 
change  on  the  gradient  or  groundwater  velocity.  /  is  estimated  as  follows: 

/  =  q'/H 
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where  q'  =  net  infiltration,  m  /  year 

H  =  leachate  penetration  depth,  m. 

The  remaining  input  parameters  and  variables  are  dependent  and  cannot  be 
generated  without  properly  matching  each  value  with  other  related  values. 
Dependencies  are  considered  with  the  objective  of  avoiding  unrealistic  or  impossible 
sets  of  data,  that  is,  to  provide  consistent  sets  of  data.  In  general,  precise  functional 
relationships  among  ail  the  dependent  variables  or  parameters  do  not  exist.  Similarly, 
observed  data  for  ail  values  taken  in  sets  do  not  exist  or  are  inadequate  in  number  to 
permit  statistical  representation  of  these  dependencies.  However,  equations  do  exist 
in  the  engineering  and  scientific  literature  to  permit  generation  of  sets  of  possible 
combinations  of  input  data.  The  parameters  and  variables  to  be  generated  as 
dependent  values  are  as  follows: 


0 


P 

V 


a 


k..k„,k. 


Ko 


dispersivities  in  the  longitudinal,  lateral  and  vertical 
directions  assumed  to  be  largely  dependent  on  distance,  x. 
porosity  of  the  soil  or  porous  media,  assumed  to  be  largely 
dependent  on  soil  p  .operties  and  parent  material, 
bulk  density  of  the  soil  or  porous  media,  largely  dependent 
on  soil  properties  including  porosity, 
groundwater  flow  velocity,  largely  dependent  on  soil 
properties  including  hydraulic  conductivities,  porosity,  bulk 
density  and  hydraulic  gradient. 

standard  deviation  of  the  gaussian  distribution 
representation  of  the  source  concentration,  related  via  mass 
balance  principles  to  leachate  volumes,  groundwater 
velocity,  porosity,  and  depth  of  leachate  penetration  into 
the  saturated  zone. 

hydrolysis  rate  constants,  dependent  on  groundwater  pH 
and  temperature,  and  on  specific  chemical  properties, 
effective  distribution  coefficient  for  each  specific  chemical. 
It  is  assumed  to  be  dependent  on  the  organic  carbon 
content  of  the  soil,  and  in  some  cases  on  the  pH  and 
specific  chemical  properties  of  the  pollutant. 


The  relationships  used  for  the  generation  of  these  parameters  and  variables  are 
discussed  below. 


g.  Dispersivity.  Guven  et  at.  (Reference  38)  reported  a  simple,  linear 
dependency  on  the  travel  distance  for  the  longitudinal  dispersivity,  of  the  form: 

Ol  =  0.093  X  +  0.007 
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where 


X 


mean  travel  distance,  m. 


Transverse  dispersivity,  a^,  has  been  studied  to  a  lesser  degree  but  its 
magnitude  is  known  to  be  less  than  the  longitudinal  dispersivity  while  maintaining 
scale  dependency  (Reference  39).  Typically,  Oj  is  related  to  by  a  simple  ratio 
leading  to  the  expression: 


Ol/Ot  =  LTR 

where  LTR  =  longitudinal/transverse  dispersivity  ratio. 

However,  such  ratios  are  often  assumed  quite  arbitrarily.  The  actual  relationship  is  the 
subject  of  current  research.  A  range  oi  LTR  has  been  reported  that  appears  to  center 
around  a  value  of  3.0,  which  has  been  selected  for  the  EPAGW  Monte  Carlo  analysis. 
For  multidirectional  flow  in  the  longitudinal  direction  the  vertical  dispersivity,  a^,  is 
quite  low.  Using  the  ratio  aja^  to  describe  vertical  dispersivity,  Gelhar  et  al. 
(Reference  40)  reported  a  range  of  30  to  1860,  with  an  average  of  400.  Because  of 
the  uncertainty  surrounding  proper  specification  of  values  for  vertical  dispersivity,  it 
is  varied  uniformly  from  40  to  400  in  the  Monte  Carlo  routine. 

The  data  generation  approach  for  dispersivity  can  be  summarized  by  the 
following  equations: 

Ol  =  0.093  X  +  0.007 

Ot  =  0.0333  X 

Oy  =  0.0025  X  -  0.01  X 

where  X  is  the  downgradient  exposure  point  distance  selected  for  the  implementation 
of  the  decision  rule. 

h.  porosity,  0.  It  is  the  ratio  of  the  volume  of  voids  of  a  given  soil  to  the 
total  volume  of  soil.  It  is  largely  a  function  of  particle  size.  For  small  particle  size  like 
clay,  porosity  increases  to  a  maximum  of  about  0.5.  Porosities  of  coarse  media  (e.g. 
gravel)  decrease  to  a  minimum  of  about  0.3.  These  measured  ranges  suggest  a  strong 
correlation  with  mean  particle  diameter,  d.  Data  reported  by  Davis  (Reference  41 ) 
were  used  to  develop  a  regression  equation  relating  porosity  to  mean  particle  size  as 


follows: 

0 

0.261  -  0.0385  In  (d ) 

where 

In  (d )  = 

natural  logarithm  of  the  mean  particle  diameter. 
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The  distribution  for  porosity  is  generated  from  a  seed  distribution  for  particle 
size  diameter.  After  a  uniform  and  a  log-uniform  distribution  for  the  particle  size 
diameter  were  investigated,  the  log-uniform  distribution  was  selected  because  it  more 
heavily  weights  the  influence  of  smaller  particle  sizes  and  because  the  related  velocity 
distribution  is  more  consistent  with  observed  data. 

i.  Bulk  density  is  defined  as  the  mass  of  dry  soil  divided  by  its  total  or  bulk 
volume.  According  to  Freeze  and  Cherry  (Reference  42): 

0  =  1  - 

Pp 

where  Pp  =  particle  size  density,  g/cm® ; 

Pb  =  bulk  density,  g/cm^ 

The  particle  density  of  soil  materials  varies  over  a  very  narrow  range,  with  an 
average  2.65  g/cm^.  Substituting  this  value  in  the  porosity  equation,  the  bulk  density 
may  be  calculated  as  a  function  of  porosity  as  follows: 


p^  =  2.  65  (  1  -  e  ) 


j.  Velocity,  V.  The  velocity  of  groundwater  is  a  major  determinant  of  the 
transport  of  solutes  in  subsurface  systems.  In  uniform  porous  media  it  is  the  dominant 
factor  and  must  be  properly  specified  in  the  Monte  Carlo  process.  Dependencies 
among  its  input  data  must  be  preserved  while  generating  realistic  values  of  velocity. 

Velocities  are  related  to  soil  properties  and  other  site-specific  factors  through 
Darcy's  Law.  Assuming  steady  flow  in  uniform,  saturated  media,  Darcy's  Law  states 
(in  simplified  form): 


where  K,  =  saturated  hydraulic  conductivity,  cm/sec 
S  =  hydraulic  gradient. 

The  saturated  hydraulic  conductivity  is  a  measure  of  the  ease  with  which  water 
flows  though  porous  media.  For  any  given  fluid,  it  is  a  function  of  the  medium 
properties  including  particle  size,  grain  shape,  connectivity  and  tortuosity.  Several 
approximate  functional  relationships  to  estimate  the  value  of  K,  have  been  proposed. 
The  most  notable  among  them  is  the  Karmen-Cozeny  equation  (Reference  43): 


47 


where  6  -  porosity; 

d  -  mean  particle  diameter. 

The  hydraulic  gradient,  S,  is  a  function  of  the  local  topography,  groundwater 
recharge  volumes  and  locations,  and  the  influence  of  withdrawals,  and  it  is  indirectly 
related  to  porous  media  properties.  Since  there  is  no  functional  relationship  to  express 
these  dependencies,  another  independent  seed  distribution  is  required  to  generate  this 
variable.  To  prevent  unrealistic  conditions  due  to  very  large  values  of  the  velocity 
caused  by  high  values  of  both  and  S,  the  velocity  field  may  L  j  bounded  such  that 
a  fixed  maximum  is  not  exceeded. 

k.  Standard  deviation  of  the  Gaussian  distribution  for  the  source 
concentration,  a.  This  parameter  reflects  the  nature  and  extent  of  the  leachate 
interaction  with  the  groundwater.  From  mass  balance  principles  it  may  be  stated  that 


a 


y/2n  Ve  HC^ 


where 

9 

= 

unit  areal  flux  of  leachate  through  the  land  disposal  facility, 

K 

m  /  year; 

area  of  disposal  facility,  m 

V 

= 

groundwater  velocity,  m  /year; 

e 

= 

saturated  zone  porosity; 

H 

= 

leachate  penetration  into  the  saturated  zone,  m; 

c. 

= 

contaminant  concentration  in  the  leachate; 

Co 

= 

contaminant  concentration  in  the  mixing  zone  beneath  the 

Assuming  that 

facility. 

the  leachate  concentration  is  the  same  as  the  maximum 

concentration  of  the  Gaussian  concentration  distribution,  a  can  be  calculated  directly 
from  the  equation  above  given  the  other  variables.  This  assumption  implies  that  the 
leachate  displaces  the  groundwater  and  dilution  begins  after  advective  transport  has 
been  initiated.  Values  for  the  chemic^J  flux,  q,  and  the  area  term,  A^,  were  generated 
by  independent  seed  distributions.  For  mathematical  reasons  (boundary  effects)  the 
constraint  that  the  ratio  H/B  be  less  than  0.5  must  also  be  made.  The  minimum 
saturated  thickness  is  set  to  3  meters. 
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I.  Hydrolysis  rates.  Hydrolysis  rates  depend  on  the  chemical  nature  of  the 
pollutant  and  have  to  be  taken  from  the  literature  or  measured  experimentally  in  the 
laboratory.  Acid  and  base-catalyzed  hydrolysis  rates  depend  on  groundwater  pH,  and 
acid,  base  and  neutral  hydrolysis  rates  are  a  function  of  temperature.  The  temperature 
dependency  has  been  described  using  the  Arrhenius  equation.  Using  the  generic 
activation  energy  recommended  by  Wolfe  (Reference  44)  of  20  kcal/mole,  the 
temperature  correction  factor  can  be  written  as  : 


=  exp 

where 

-  second-order  hydrolysis  rate  constant  for  acid,  neutral,  or 

base  conditions  at  temperature  T; 

=  second-order  hydrolysis  rate  constant  for  acid  neutral  or 

base  conditions  at  reference  temperature  T,) 

T 

=  temperature, 

T, 

=  reference  temperature,  °K. 

m.  Distribution  Coefficient.  In  most  cases,  the  sorption  process  is  dominated 
by  hydrophobic  binding,  it  is  possible  then  to  reiate  the  distribution  coefficient  directly 
to  the  soil  organic  carbon.  As  stated  earlier,  the  dependency  is  given  by: 

ko  =  (koc)  ( 


where 

koc  =  distribution  coefficient  normalized  to  organic  carbon; 
f^  =  fractional  organic  carbon. 

The  values  for  fractional  organic  carbon  are  generated  as  an  independent 
parameter.  For  other  binding  mechanisms,  as  those  presented  by  polar  or  ionizable 
compounds,  adjustments  are  made  on  a  case-by-case  basis  (Reference  45).  In  cases 
where  reliable  relationships  do  not  exist,  measurements  are  required. 


8.  EPASF,  Modified  ERA  Model  For  Analysis  Of  Impact  On  Surface  Water 

As  part  of  rules  for  land  disposal  restrictions  proposed  by  the  U.S. 
Environmental  Protection  Agency  (Reference  20),  models  were  developed  for  the  back- 
calculation  of  screening  levels  of  constituents  in  land  disposal  systems  under 
conditions  of  uncertainty.  The  models  developed  for  this  purpose  considered  the 
impacts  of  hazardous  constituents  both  in  groundwater  and  in  surface  waters 
contaminated  via  migration  through  groundwater.  These  models  have  been  modified 
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for  down-gradient  calculation  of  concentrations  under  uncertainty  and  have  been 
included  in  the  workstation  Advisory  System  under  the  acronyms  EPAGW  (previous 
section)  and  EPASF.  EPASF  addresses  the  contamination  scenarios  concerned  with 
impacts  in  and  through  surface  waters.  These  are: 

Scenario  2A:  Violation  of  an  environmental  standard  in  surface  water 
contaminated  by  leachate.  This  scenario  involves  failure  of  the  waste 
container,  followed  by  transport  in  groundwater  and  mixing  with  uncontaminated 
stream  water. 

Scenario  2B:  Exposure  of  humans  through  drinking  water  coming  from 
surface  water  contaminated  by  leachate.  This  continues  the  previous  scenario  with 
downstream  transport,  intake  by  water  treatment  plant,  and  exposure  of  humans  to 
the  contaminant  via  drinking  water. 

Scenario  3:  Exposure  of  humans  through  consumption  of  fish  from  surface 
waters  contaminated  by  leachate  carried  through  groundwater.  This  proceeds  as  in 
Scenario  2A,  followed  by  uptake  of  contaminant  by  fish  (bioconcentration  and/or 
biomagnification),  and  human  exposure  via  fish  consumption. 

All  three  of  these  scenarios  require  calculation  of  diluted  instream 
concentrations.  This  may  then  be  equated  to  standards  of  human  exposure  through 
drinking  water  and  human  exposure  through  consumption  of  fish,  by  appropriate 
assumptions  regarding  daily  water  intake  and  daily  average  fish  consumption.  For 
instance,  the  pathway  analysis  through  fish  consumption  is  based  on  an  average 
consumption  rate  of  6.5  grams  of  fish  per  day. 

To  derive  an  approximate  solution  we  first  assume  that  the  average 
concentration  at  the  ground  water  outlet  to  the  surface  water,  C,,  can  be  related 
linearly  to  the  leachate  concentration,  C^,  as: 


where  z,  is  a  groundwater  attenuation  factor. 

It  is  assumed  that  the  saturated  zone  transport  has  reached  steady-state.  The 
contaminant  mass  flux  leaching  from  the  site  into  the  groundwater  system,  is 
given  by: 

where  is  the  volumetric  rate  of  leaching  from  the  site,  and  Q  is  the  leachate 
concentration.  Solute  transport  in  groundwater  from  the  site  to  the  stream  results  in 
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a  plume  that  intercepts  the  stream  over  an  area  with  an  average  concentration  C,. 
C,  corresponds  to  the  average  of  the  actual  point  concentration,  which  is  assumed  to 
be  a  Gaussian  distribution  over  the  effective  flow  area.  If  the  groundwater  steady- 
state  seepage  velocity  is  represented  by  V^,  the  contaminant  mass  flux  from  the 
groundwater  system  into  the  surface  water  system,  is  given  by: 

"i ' 

where  Q,  is  the  contaminated  ground  water  discharge  rate. 

At  steady  state,  m,  is  related  to  /nj,  by: 


ny  = 


where  is  the  fraction  of  the  contaminant  mass  not  transformed  by  hydrolysis  or 
initial  speciation  in  the  groundwater  (as,  at  steady  state,  lateral  dispersion  of  the 
plume  does  not  affect  the  total  mass  loading  to  the  stream). 


Combining  these  equations  gives  an  expression  for  the  concentration  dilution 
factor  due  to  transport  in  groundwater: 


Q  = 


The  parameters  Q,,  and  are  estimated  in  the  simplest  case  from: 

_  ^(1 

^  ■  (86400)  (365.25) 


fn  =  exp  {-/^^) 

where: 

Ql  =  rate  of  percolation  through  the  land  disposal  unit,  m^/sec. 
P  =  average  annual  precipitation  rate,  mlyear. 

f^^  =  runoff  fraction. 

=  surface  area  of  the  waste  site.  m^. 

/C,  as  total  effective  decay  constant  in  groundwater,  1/yrs. 

r,  =  time  taken  by  the  contaminant  to  travel  from  the  land 

disposal  unit  to  the  stream  entry  point,  years. 


51 


The  factor  is  estimated  exactly  as  in  the  model  EPAGW.  Similarly,  the  travel 
time  of  the  constituents  in  groundwater  is  given,  in  the  one-dimensional  case  without 
dispersion,  by: 


T 


9  '  Dg 


where: 


distance  from  site  to  stream,  m. 

groundwater  seepage  velocity,  m  /  yr. 

fraction  of  compound  that  is  dissolved,  calculated  as  in  the 

model  EPAGW. 


When  the  contaminated  ground  water  enters  the  stream  it  mixes  with  surface 
water  supplied  by  the  upstream  watershed.  Assuming  that  the  downstream  end  of 
the  plume  entry  is  given  by  x«0,  and  that  lateral  concentration  gradients  disappear 
due  to  mixing,  the  laterally  averaged  concentration,  C.  increases  with  x  reaching  a 
maximum  near  the  downstream  end  of  the  impacting  plume  (xsQ). 

From  mass  balance  considerations,  at  x^O, 

,aL 


C. 


where  Q*,  is  the  stream  base  flow  at  x=0,  given  by: 

*  "  (86400)  (365.  25) 


where  A,  is  the  surface  area  of  the  upstream  watershed  and  it  is  assumed  that  the 
average  annual  precipitation  rate,  P.  and  the  average  runoff  coefficient,  f^,  are  the 
same  for  the  waste  site  and  the  entire  watershed.  Combining  equations  then  yields 
the  laterally  averaged  concentration  at  the  downstream  edge  of  the  groundwater 
plume,  C,,  as: 

r*  _  ^  H^wr* 

- A — 


The  steady-state  laterally  averaged  value  of  concentration  downstream  may 
then  be  approximated  through  an  attenuation  factor,  Zy  =  O'*,  where  &  =  KxAJ,  K  - 
decay  rate  constant,  sec'\  and  U  =  mean  downstream  velocity,  m/sec.  In-stream 
decay  processes  include  sorption,  hydrolysis  and  volatilization.  Hydrolysis  is 
calculated  dependent  on  pH  and  organic  carbon  fraction  of  the  suspended  sediment, 
in  a  manner  similar  to  that  employed  for  hydrolysis  in  EPAGW. 


As  programmed,  the  model  provides  considerable  flexibility  through  the 
determination  of  f^,  the  mass  loading  factor,  and  r,,  the  time  taken  by  the  contaminant 
to  travel  from  the  land  disposal  unit  to  the  stream  entry  point.  The  situation  described 
above  demonstrates  the  one-dimensional  flow  case  without  dispersion.  Options  are 
also  included  to  derive  and  r,  with  the  inclusion  of  dispersion,  using  a  one- 
dimensional  advection-dispersion  solution,  and  by  use  of  a  three-dimensional  transport 
equation.  For  the  latter  case,  the  groundwater  transport  equation  used  in 
determination  of  is  the  same  as  is  documented  for  the  model  EPAGW. 


9.  LTIRD,  Radial  Flow  Disposal  Systems 

The  advection-dispersion  equation  for  plane  radial  flow  (Reference  43) 
n*ay  be  written  as: 


1  a  .  ac 

r  dr  [  dr\  dr  dt 


For  steady  plane  radial  flow  (but  transient  mass  transport),  replacing  the  dispersion 
coefficient  D  by  0|,V,  the  following  expression  is  obtained  (Reference  16): 


Q,V^  -  ^  ^ 

^-2  dr  dt 


dr 


Consider  a  confined  aquifer  with  thickness  b  being  recharged  through  a  fully 
penetrating  well  at  a  constant  rate  Q.  If  the  concentration  of  a  chemical  in  the 
recharge  fluid  is  Cq  (and  the  concentration  of  that  chemical  in  the  aquifer  water  was 
originally  zero),  the  equation  above  may  be  rewritten  as: 


1  ^  I  dCp  _  dCp 

^  D  dRo  ^  D  D  D 


where: 


/  = 

^  2nbnctL 
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Initial  and  boundary  conditions  for  this  problem  are: 


^  f  ^ 

lim  Cj,  (r^,  tj,)  =  0 

ro-^m 

where  is  the  dimensionless  well  (source)  diameter.  The  solution  in  the  Laplace 
transform  domain,  in  terms  of  Airy  functions,  is: 


where: 

Z-D  =  Laplace  transform  of  dimensionless  concentration, 

a  =  Laplace  transform  parameter. 

Y  =  {sfo  +  %) 

Yo  =  (sro^  +  %) 

The  general  form  of  the  Airy  function  can  be  found  In  Abramowitz  and  Stegun 
(Reference  46).  An  asymptotic  expansion  yields: 

Ai  (z)  •  ^  y/n  z'^^*  5^  (“1)*  5'* 

for  |z|  <  1 

where: 


The  computer  code  LTIRD  calculates  the  dimensionless  concentration  of  a 
solute  injected  into  an  aquifer,  it  has  been  adapted  for  single  rotary  distributors  in  the 
workstation  Advisory  System,  but  should  not  be  used  for  multiple  distributors. 
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10.  RESSQ,  Semi-analyticai  Model  Based  On  Complex  Velocity  Potential 

Semi-analytical  methods  are  more  powerful  than  analytic  methods  in 
terms  of  representation  of  the  flow  regime,  and  simpler  than  most  of  the  complete 
numerical  techniques.  These  methods  apply  a  well-known  concept  of  fluid  mechanics: 
the  complex  velocity  potential.  A  major  limitation  is  that  they  apply  only  to 
steady-state  two-dimensional  fluid  flow  through  homogeneous  media.  The  computer 
program  RESSQ  calculates  two-dimensional  contaminant  transport  by  advection  and 
adsorption  (no  dispersion  or  diffusion)  in  a  homogeneous,  isotropic  confined  aquifer 
of  uniform  thickness  when  regional  flow,  sources  and  sinks  create  a  steady-state  flow 
field  (Reference  16). 

The  velocity  potential  is  generally  defined  as: 

=  Kh  *  c 


Therefore,  a  component  of  the  specific  discharge  or  Darcy  velocity  vector  in  any 
arbitrary  direction  x  is: 


Qjc 


The  stream  function  can  be  obtained  with  a  known  velocity  potential  by  using 
the  Cauchy-Riemann  equations: 

60  _  dijr 
dx  dy 


60  ^  _  6\|f 

By  dx 


Both  the  stream  function  and  velocity  potential  are  harmonic  functions  in  that 
they  satisfy  the  Laplace  equation:  therefore,  their  application  is  restricted  to 
steady-state  planar  flow  fields.  The  complex  velocity  potential  of  a  uniform  flow  with 
Darcy  velocity  U,  in  a  direction  making  an  angle  a  with  the  positive  x  axis,  is: 

W=  -UZe-‘°  +  c 
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Substituting  for  complex  numbers  Z  and  the  velocity  potential  and  the 
stream  function  for  such  a  flow  system  are  obtained: 

ly*  0  +  /  ijr  =  -4^  jr  +  ^ )  ( cos  a  -  /  sin  a)  *  c 
0  =  -O(xcosa  ♦  y  sin  a)  =  constant 
0  =  4>(jirsina  ~  y  cos  o)  =  constant 


The  latter  two  equations  represent  equipotentials  and  streamlines,  respectively. 
Components  of  specific  discharge  are: 


<7,  -  =  Ueos  a 


The  complex  velocity  potential  of  a  source  with  strength  m  located  at  the  point 

Zois: 


m  \r\  {Z  -  Zq)  *  c 


If  the  source  represents  a  well  which  is  being  recharged  at  rate  Q  into  an  aquifer  of 
thickness  b,  then  the  strength  m  of  the  source  is  simply: 


/w  = 


Q 

2nb 


Substituting  for  complex  numbers  Z  and  Zo.  the  velocity  potential  and  stream 
function  are: 


-Q 


In  -  / 


7  _  /  O 
2nb 


tan 


-1 


y  -  /o 


X  -  Xa 


+  C 


0  =  -^ln  [{X  -  xy  *  (y  -  y^)^]*  c, 


0  = 


-  ->ro) 


+  C2 


where  Xg  and  /o  sre  the  coordinates  of  the  source  and  x  and  y  are  the  coordinates  of 
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a  point  where  the  velocity  potential  and  stream  function  are  calculated.  Components 
of  specific  discharge  based  on  the  above  definitions  are: 

Sx  2S5  (x-Ai)*  .(/-yo)* 


(r  -  yj 


By  Snb  (jr  -  Jfo)®  *(/ -  yo)* 


The  complex  velocity  potential  for  a  positive  doublet  is: 

w=  -2^  +  c 
|2|» 


where 

Z  =  conjugate  of  complex  number  Z  and 
\Z\  =  modulus  of  complex  number  Z 


Substituting  for  these  two  values,  the  velocity  potential  and  stream  function  for 
a  positive  doublet  located  at  the  origin  are: 

x^  * 

0  =  -  A  =  constant 

X2  +  y2 

t|r  =  =  constant 

^  x^  * 


The  latter  two  equations  can  be  rearranged  to  give: 


The  first  equation  represents  equipotentials  and  describes  circles  with  centers  along 
the  X  axis.  Similariy,  the  second  equation  represents  a  group  of  circles  with  centers 
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on  the  y  axis  and  tangent  to  the  x  axis.  These  circles  are  the  streamlines  for  such  a 
doublet. 

Since  the  Laplace  equation  is  a  linear  partial  differential  equation,  the  principle 
of  superposition  applies:  as  many  flow  components  as  needed  can  be  superimposed 
to  obtain  the  expression  for  the  complex  velocity  potential  of  an  entire  system.  For 
example,  one  or  several  point  sources  of  contaminant  recharge,  together  with  some 
groundwater  discharge  wells,  can  be  combined  with  a  uniform  regional  groundwater 
flow  regime.  The  overall  complex  velocity  potential  may  be  written  as: 

where: 

W  s  overall  complex  velocity  potential  of  the  system. 

U  =  Oarcy  velocity  of  uniform  regional  flow. 

a  «  angle  between  the  direction  of  regional  flow  and  the  positive  x 

axis. 

b  «  aquifer  thickness. 

Qy  =  rate  of  discharge  from  well  /. 

=  rate  of  discharge  from  well  k. 

The  velocity  potential  of  the  above  system,  the  real  part  of  W,  is  (Reference 

16): 

=  -i/(  JT  cosa  *  /  sin  a)  *  T  In  *  (/-//)*] 


and  the  expression  for  the  stream  function,  the  imaginary  part  of  W,  becomes: 
t|r  =  6(  jr  sin  a  -  /  cosa)  +  T  tan  ’ 


y-Vk 


At  any  given  point  with  coordinate  (x,y),  components  of  the 
specific  discharge  for  the  overall  system  may  be  written  as: 
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Qx 


6/ cos  a  - 


_ ( x~Xj) 

h  (x-Xj)^  *  {y-yj)^ 


h  2wZ>  (jr-jf*)2  +  (x-yj* 


Qy  = 


=  -^  = 


ay 


4/sin  a  - 


/V 


( /-//) 


2JI/>  (Jf-Jfy)*  H 

4  ( y-yi>) 

2irA  (jf-jr^)*  + 


( y-y/) 


Components  of  the  average  pore  water  velocity  for  an  individual  fluid  particle 
moving  through  the  overall  flow  system  are  (by  introducing  a  retardation  factor,  /?): 


^cx  = 


Si. 

nR 


cy 


Si 

nR 


The  path  line  travelled  by  a  contaminant  particle  can  be  divided  into  increments 
dl,  traversed  in  time  intervals  dt.  The  projections  of  d!  on  the  x  and  y  axes  are  given 
by  dx  and  dy,  respectively: 


dx  =  v^dt  = 


Qxdt 

nR 


dy  =  v^dt  = 


dydt 

nR 


d!  -  V(  dx^  *  dy^)  =  ,j(  <7/  *  <7/)  ^ 


Numerical  integration  of  the  equation  for  d!  yields  travel  time  between  any  two 
points  of  a  given  streamline.  If  a  contaminant  particle  is  at  a  point  (Xy,  y^  at  time  t,  its 
new  position  at  time  f  +  r/f  on  the  same  streamline  can  be  calculated  by  using: 
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According  to  Javandel  etaL  (Reference  1 6),  the  combination  of  a  uniform  flow 
in  the  positive  x  direction  with  a  positive  doublet  and  a  point  source  (both  centered 
at  the  origin)  represents  outflow  from  a  completely  penetrating  cylindrical  pond  in  the 
presence  of  a  uniform  flow  in  the  positive  x  direction.  The  complex  velocity  potential, 
then,  is: 


W=  -UZ* 


QZ 

1-21" 


i.ln 


2nl} 


Z*  c 


and  the  velocity  potential  and  stream  function  are: 


0  =  -Ux  + 


Ox 


4 


Attb 


In  (x*  +  y^)  + 


ijr  =  ~Uy  + 


Sk—  -  _S. 


JK-2  ^  y2 


2nb 


tan 


*  Co 


where: 

U  =  Darcy  velocity  of  uniform  flow  in  the  positive  x  direction. 

Qp  =  rate  of  outflow  from  the  pond. 

b  =  thickness  of  the  aquifer. 

Q  s:  constant  of  the  doublet. 


The  value  of  the  constants  in  the  velocity  potential  equation  can  be  determined 
such  that  the  velocity  potential  satisfies  the  appropriate  boundary  conditions.  Holding 
0  constant  and  equal  to  HqZXt  =  /q. 


0  =  -  4/x  + 


t/yjr 


Q> 


Anb 


In 


x^  +  y^ 


Incorporating  the  velocity  potential  of  sources  and  sinks,  the  result  is: 


n-  Ux* 


A*" 

Au  b 


-  yj^ 


S 


4 

Anb 


In 


(jr-jr^)‘ 


^  ky-yj^^ 

*  yk 


where  Qy  and  Q,,  are  the  rates  of  discharge  and  recharge  of  sinks  and  sources, 
respectively.  Components  of  the  average  pore  water  velocity  at  any  point  {x.y)  within 
the  overall  flow  regime  where  the  velocity  potential  is  defined  may  be  written  as: 


_  1  _  4/  .  ^^0 


n  dx 


x^ 


A. 


2nnb  +  y^ 


(j^2  +  y2j2 

^  __g _ ^  ^  _ jJ^-J^k) 

H  {x-xj)^  +  (y-y/)*  hi  27r/7Z>  {x-x^^  +  {y-yy 


V  =  -1  = 

^  n  dy  n 


and 

Zxy 


Op 


N 


y'  ^ _ M  ip _  ^  ^ _ V  y  Xk! _ 

h  {x-xj)^  +  (y-yy)^  M  271/7/7  (jr-jr^)2  +  (y-y^)^ 


(y-yy) 


(j^2  +  y2j2 

"  4 


Znnb  x^  + 

(y-y-r) 


11.  MOC,  Method  of  Characteristics  Two-Dimensional  Solute  Transport 

Model 

This  model  can  be  applied  to  a  wide  variety  of  one-  and  two-dimensional 
problems  involving  steady-state  or  transient  flow.  It  computes  changes  in 
concentration  due  to  the  following  processes  (References  21,  22  and  23)  : 

a.  convective  transport  by  which  dissolved  chemicals  move  with  the 
flowing  groundwater; 

b.  hydrodynamic  dispersion,  by  which  molecular  and  ionic  diffusion 
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and  smali'Scale  variations  in  the  velocity  of  flow  through  the  porous  media  cause  the 
path  of  dissolved  molecules  or  ions  to  spread  from  the  average  direction  of 
groundwater  flow; 


c.  fluid  sources,  causing  mixing  or  dilution; 

d.  reactions,  by  which  the  concentration  of  chemical  is  modified  by 
its  interaction  with  other  species  present  in  the  groundwater  solution  or  the  solid 
aquifer  materials. 

The  model  assumes  that  the  gradients  of  fluid  density,  viscosity  and 
temperature  do  not  affect  the  velocity  distribution.  The  solute  can  be  reactive  or 
conservative,  and  the  aquifer  can  be  heterogeneous  and/or  anisotropic.  The  computer 
program  first  solves  the  governing  equation  that  describes  the  transient  twu- 
dimensional  areal  flow  of  a  homogeneous  compressible  fluid  through  a 
nonhomogeneous  anisotropic  aquifer: 


a/yy 

dx/ 


w 


where 

Tf  =  transmissivity  tensor  (IVT)  =  KJD 

Kf  -  hydraulic  conductivity  tensor  {L/D 

S  =  storage  coefficient 

W  =  source  or  sin',  term 

Xi  =  cartesian  coordinates  {L) 

h  =  hydraulic  head  {L) 

b  =  aquifer  thickness  (I). 

The  source  or  sink  term,  W,  is  the  volume  flux  per  unit  area.  When  only  the 
following  terms  are  considered:  (a)  direct  withdrawal  or  pumpage  (well  pumpage, 
evapotranspiration  or  well  injection);  (b)  steady  leakage  into  or  out  of  the  aquifer 
through  a  confining  layer,  then  the  term  W  may  be  expressed  as 


=  0(x.y,t)  - 


where 

Q  =  rate  of  withdrawal  (positive  sign)  or  recharge  (negative  sign), 

(L/7); 

Kg  =  vertical  hydraulic  conductivity  of  the  confining  layer,  (L/7);  m  - 

thickness  of  the  confining  layer,  (Z.); 
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H,  «  hydraulic  head  in  the  source,  (Z.). 


Darcy's  Law  provides  a  basis  to  calculate  the  average  seepage  velocity  of 
groundwater: 


^  dh 

€  dXj 


where 


V,  s  seepage  velocity  in  the  direction  of  x,  (t/7); 

€  »  effective  porosity  of  the  aquifer. 

The  mass  transport  equation  is: 


where 

C  =  concentration  of  soiute  {M/L^); 

Dg  =  coefficient  of  dispersion  (Z.V7); 

=  soiute  concentration  in  source  or  sink  fluid  [M/L\ 

Rg  =  rate  of  addition  or  removal  of  soiute  by  physical  or  chemical 

reactions  {M/L^D; 

b  =  saturated  thickness  of  the  aquifer  (L). 

The  first  term  on  the  right-hand  side  of  the  transport  equation  describes  the 
hydrodynamic  dispersion,  the  second  term  represents  the  convective  transport,  and 
the  third  and  fourth  term  describe,  respectively,  fluid  sources  or  sinks  and  changes  in 
concentration  due  to  chemical  or  physical  reactions  occurring  in  the  groundwater 
solution.  The  first  step  in  solving  this  equation  is  to  estimate  the  value  of  Dg. 
According  to  Scheidegger  (Reference  47),  the  dispersion  coefficient  can  be  computed 
as  a  function  of  the  velocity  of  groundwater  flow  and  the  nature  of  the  aquifer. 
Assuming  that  the  molecular  diffusion  is  negligible. 


D  -a 

and 

\y\  ♦  If 
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v^tare 


dispersivity  of  the  aquifer  {L); 

components  of  velocity  in  the  m  and  n  directions, 
respectively  (4/7). 


For  an  isotropic  medium,  the  dispersivity  tensor  may  be  defined  in  terms  of  the 
longitudinal  and  transverse  dispersivities  of  the  aquifer  (o^  and  a^,  respectively).  The 
components  of  the  dispersion  coefficient  may  be  described  as: 

^2  y2 

^  W  *  “’’M 

y2  y2 

I  y^2 


where  Dl  -  V\  and  ~  Or  I  are  the  longitudinal  and  transverse  dispersion 
coefficients. 

The  sequence  of  steps  involved  in  solving  these  equations  is  to  determine  the 
velocity  distribution  from  the  flow  equation  and  Darcy's  Law.  Next,  given  the  values 
of  Oi  and  a^,  Darcy's  Law  is  solved  subject  to  certain  boundary  and  initial  conditions. 

Flow  and  mass  transport  equations  are  treated  as  uncoupled  equations  in  that 
concentration  changes  do  not  affect  the  flow  field;  this  applies  where  density 
differences  are  negligible,  which  is  the  case  in  most  contaminant  problems,  and  one 
of  the  assumptions  in  the  MOC  model.  There  are  two  situations  where  solutions  to 
the  previous  equations  may  be  applied:  (1)  to  assess  the  impact  of  proposed 
subsurface  waste  disposal  sites  that  have  not  yet  been  contaminated;  (2)  to  assess 
the  impact  of  already  contaminated  sites,  predict  plume  migration  and  recommend 
remedial  actions. 

Numerical  Solution 

Exact  analytical  solutions  to  the  partial  differential  equations  of  flow  and  solute 
mass  transport  cannot  be  obtained  directly  due  to  vari:  properties  of  aquifers  and 

complex  boundary  conditions.  The  solutions  must  be  oximated  by  a  numerical 
scheme.  The  basic  method  is  to  break  up  continuous  space  into  cells,  approximate 
the  governing  partial  differential  equations  by  differences  between  the  values  of  the 
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parameters  over  the  network  and  then  compute  the  discrete  parcel.  MOC  utilizes  a 
rectangular,  uniformly  spaced,  block-centered,  finite-difference  grid— in  which  nodes 
are  defined  at  the  centers  of  the  rectangular  cells. 

The  numerical  solution  technique  first  solves  the  equation  describing  the 
transient  two-dimensional  flow  of  a  homogeneous  fluid  through  a  nonhomogeneous 
anisotropic  aquifer  with  the  following  implicit  finite-difference  approximation: 


r  .  T  ^l*\.J.k-^l,J.k 

^  r  ^IJ-yk-^*l,l.k  .  <r  ^IJ*yk-^l.j.k 

=  5  ^i.J.f^i.J.k-\  ^  ^  tu  ^  ^  \ 

(A/)  LxLy 

where  i,j,k  =  indices  in  the  x,y,  and  t  dimensions,  respectively; 

Ax,  Ay,  Af  =  increments  in  the  x,y  and  time  dimensions,  respectively; 

-  volumetric  rate  of  withdrawal  or  recharge  at  the  (4/)  node  (Z.V7). 


This  equation  is  solved  numerically  for  each  node  in  the  grid  using  an  iterative 
alternating-direction-implicit  (ADI)  procedure  (References  24,  48,  and  49).  After  the 
head  distribution  is  computed,  the  velocity  of  the  groundwater  flow  is  computed  for 
each  node  using  an  explicit  finite  ui'^erence  form  of  Darcy's  equation: 


2bx 


Next,  the  solute  transport  equation  is  solved.  This  equation  describes  the  two- 
dimensional  areal  transport  and  dispersion  of  a  given  reactive  dissolved  chemical 
species  in  flowing  groundwater.  This  equation  is  solved  using  the  method  of 
characteristics,  through  a  three  step  procedure.  If  saturated  thickness  is  considered 
as  a  variable  and  the  convective  transport  term  is  expanded,  the  mass  transport 
equation  may  be  rewritten  as 


dC  1  a 
dt  b  dx, 


dC 

dXj 


1/  J£  .  ]/  Jl . 

'  dx,  y  dXf 


Rk 


which  is  in  the  form  solved  by  the  computer  program. 


Changes  with  time  in  properties  of  the  fluid  (e.g.,  concentration)  may  be 
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described  either  for  fixed  points  within  a  stationary  coordinate  system  as  succesive 
particles  pass  the  reference  points,  or  for  reference  fluid  particles  as  they  move  along 
their  respective  paths  past  fixed  points  in  space.  The  derivative  dC/df  is  the  rate  of 
change  of  concentration  as  observed  from  a  fixed  point,  whereas  dC/dt  is  the  rate  of 
change  as  observed  when  moving  with  the  fluid  particle  (material  derivative). 

e  Particle  tracking.  This  first  step  solved  for  the  change  in  concentration 
over  distance.  It  consists  of  placing  a  number  of  traceable  particles  in  each  ceil  of  the 
finite-difference  grid,  to  form  a  set  of  points  that  are  distributed  in  a  geometrically 
uniform  pattern  throughout  the  area  of  interest.  The  location  of  each  particle  is 
specified  by  its  x  and  y  coordinates.  The  initial  concentration  assigned  to  each  point 
is  the  initial  concentration  at  the  node  of  the  cell  containing  the  particle.  For  each  time 
step  every  point  is  moved  a  distance  proportional  to  the  length  of  the  time  increment 
ind  the  velocity  at  the  location  of  the  point.  The  new  position  of  the  particle  is  then 
computed  and,  after  all  points  have  been  moved,  the  concentration  at  each  node  is 
temporarily  assigned  the  average  of  the  concentrations  of  all  points  located  within  the 
area  of  that  cell. 

A  two-step  explicit  finite-difference  approximation  is  now  used  to  calculate  the 
new  nodal  concentrations  at  the  end  of  the  time  period.  The  changes  in  concentration 
caused  by  hydrodynamic  dispersion,  fluid  sources,  divergence  in  velocity,  changes  in 
saturated  thickness,  adsorption  or  chemical  reaction  are  calculated  using  an  explicit 
finite-difference  approximation.  This  change  in  concentration  can  be  considered  as 
the  sum  of  two  separate  terms: 

^Qj.k  ~  *  i.^Q,/,k)k 


where  the  subscript  /  represents  the  change  in  concentration  caused  by  hydrodynamic 
dispersion,  and  the  subscript  //  is  the  change  caused  by  an  external  fluid  source, 
changes  in  saturated  thickness,  adsorption  or  chemical  reaction: 


i^Q./.k)  / 


^^.i,k 

*  ^lj,k 


Later  modifications  (References  22  ar«d  23)  include  first-order  radioactive  decay 
{M  ,  and  retardation  factors  {R,)  for  linear  sorption,  non-linear  sorption  and  ion 
exchange. 

Decay:  X  = 

*  1/2 
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whert  tia  i*  the  half-life  of  the  solute,  (7).  For  example. 

Linear  Sorption  :  /^  =  1  ♦ 


Nonlinear  Sorption  :  /^(  C)  =  1  ♦  ^ 

The  decay  is  applied  directly  to  the  tracer  particles  (rather  than  at  the  nodes  of  the 
finite  difference  grid): 

dp  =  dp^  exp  (  -AA/) 

where 

=  solute  concentration  of  the  tracer  particle,  and 
k  s  time  dimension  index. 

The  exponential  decay  formulation  has  no  numerical  stability  restrictions,  but  some 
numerical  accuracy  is  lost  if  the  half-life  is  much  smaller  than  the  time  step  for  solving 
the  transport  equation. 

e  Stability  criteria.  The  explicit  numerical  solution  of  the  solute-transport 
equation  has  some  stability  criteria  associated  with  it,  which  may  require  that  the  time 
step  used  to  solve  the  flow  equation  be  subdivided  into  a  number  of  smaller  time 
increments  to  accurately  solve  the  solute-transport  equation.  The  stability  criteria  may 
be  stated  as  follows: 


Min 
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where  y  is  the  fraction  of  the  grid  dimensions  that  particles  are  allowed  to  move  at 
each  time  step  (0  i  y  i  1). 

If  the  time  step  used  to  solve  the  flow  equation  exceeds  the  smallest  of  the  time 
limits  determined  by  the  above  equations,  then  the  time  step  will  be  subdivided  into 
the  appropriate  number  of  smaller  time  increments  required  for  solving  the  solute- 
transport  equation. 

•  Boundary  and  initial  conditiona..  Several  different  types  of  boundary 
conditions  ca^'  ^  >  incorporated  into  the  solute  transport  model.  Two  general  types 
have  been  used  in  this  model:  constant  flux  and  constant-head  conditions.  A 
constant  head  boundary  can  be  used  to  represent  aquifer  underflow,  well  withdrawals 
or  well  injection.  A  finite  flux  is  defined  by  specifying  the  flux  rate  as  a  well  discharge 
or  injection  rate  for  the  appropriate  nodes.  A  no-flow  boundary  is  a  special  case  of 
a  constant  flux  boundary.  The  numerical  procedure  requires  that  the  area  of  interest 
be  surrounded  by  a  no-flow  boundary.  No-flow  boundaries  can  also  be  located 
elsewhere  in  the  grid  to  represent  natural  barriers  to  groundwater  flow.  No-flow 
boundaries  are  designated  by  setting  to  zero  the  transmissivity  at  appropriate  nodes. 

A  constant-head  boundary  represents  parts  of  the  aquifer  such  as  recharge 
boundaries  or  areas  beyond  the  influence  of  hydraulic  stress.  Constant-head 
boundaries  are  represented  by  adjusting  the  leakage  term  at  appropriate  nodes.  If  a 
constant-flux  or  constant-head  represents  a  fluid  source,  then  the  chemical 
concentration  in  the  source  fluid  must  also  be  specified.  The  initial  conditions  can  be 
determined  from  field  data  and  from  previous  simulations.  The  head  and  concentration 
in  the  aquifer  at  the  start  of  the  simulation  must  be  specified,  because  solute  transport 
depends  directly  upon  hydraulic  and  concentration  gradients. 

Monte  Carlo  Simulation  Option 

The  MOC  model  as  described  in  the  preceding  sections  involves  a  large  degree 
of  uncertainty  originated  by  the  combination  of  model  error,  natural  and  parameter 
uncertainties.  As  an  alternative  to  the  deterministic  approach  in  which  detailed  data 
are  required  for  the  simulation,  the  Monte  Carlo  technique  enables  a  prediction  which 
incorporates  the  uncertainty  associated  to  the  inputs  and  parameters  of  the  model. 
This  is  achieved  by  a  random  generation  of  the  most  sensitive  inputs  to  the  model, 
followed  by  a  large  number  of  computations  to  yield  a  well  defined  distribution  of 
outputs. 


1 2.  HELP,  Hydrologic  Evaluation  Of  Landfill  Performance 

The  Hydrologic  Evaluation  of  Landfill  Performance  (HELP)  model  was 
developed  to  facilitate  rapid,  economical  estimation  of  the  amounts  of  surface  runoff, 
subsurface  drainage,  and  leachate  that  may  be  expected  to  result  from  the  operation 
of  a  wide  variety  of  possible  landfill  designs  (Reference  51).  These  phenomena  arise 
from  the  interaction  of  a  large  number  of  complex  hydrologic  processes,  including 
precipitation,  surface  storage,  runoff,  infiltration,  percolation,  evapotranspiration,  soil 
moisture  storage  and  lateral  drainage.  HELP  takes  what  is  essentially  a  water-balance 
approach  to  the  problem,  within  a  quasi-two  dimensional  process.  The  model  uses 
climatologic,  soil,  and  design  data  to  produce  estimates  of  water  movement  across, 
into,  through,  and  out  of  landfills.  To  accomplish  this,  daily  precipitation  is  partitioned 
into  surface  storage  (snow),  runoff,  infiltration,  surface  evaporation,  evapo¬ 
transpiration,  percolation,  stored  sol!  '*^<'<sture,  and  subsurface  lateral  drainage  to 
maintain  a  water  budget.  Calculations  may  be  performed  for  up  to  nine  layers  of  a 
landfill  design.  These  layers  may  include  a  vegetative  layer,  other  types  of  vertical 
percolation  layers,  lateral  drainage  layers,  barrier  soil  layers  and  waste  layers. 

Surface  Runoff.  Surface  runoff  is  computed  by  the  SCS  curve  number 
technique  (Reference  51).  This  method  was  chosen  by  the  authors  because  it  is  well 
established,  easy  to  use,  and  the  required  input  data  is  generally  available.  Generally, 
the  curve  number  for  a  watershed  is  determined  for  average  moisture  conditions  in  the 
SCS  method  (CA/J.  The  curve  number  for  the  lowest  antecedent  moisture  conditions, 
CN^  is  related  to  CN„  by  a  polynomial  developed  for  the  CREAMS  model  (Reference 
52): 

CN,  =  -16.  91  +  1.  348  CN„  -  0.  01379  CNy  2  +  0.  0001177  CN„  ® 


The  maximum  retention  factor  for  a  soil,  S^ 

<5  _  1000 
^  ■  "CFIT 


is  then  determined  from: 
-  10 


From  this  information  we  can  calculate  the  daily  depth  weighted  retention  factor,  Sj, 
and  the  daiiy  runoff,  Q,-,  using  the  method  documented  by  Knisel  (Reference  52).  The 
soil  profiie  of  the  vegetative  or  evaporative  depth  was  divided  into  seven  segments. 
The  thickness  of  the  top  segment  was  set  equai  to  1  /36  of  the  thickness  of  this  depth, 
whiie  the  thickness  of  the  second  segment  was  5/36  of  this  depth.  Each  of  the 
remaining  5  segments  was  defined  as  1/6  of  the  thickness  of  the  vegetative  or 
evaporative  depth.  We  may  then  state: 


S,  = 
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where 


soil  water  content  of  segment  j,  inches 
saturated  capacity  of  segment  j,  inches 
wilting  point  of  segment  j,  inches 


where: 

Dj  a  depth  to  bottom  of  segment  j,  inches 
VD  a  vegetative  or  evaporative  depth,  inches, 
and 

Q  = 


where: 

P  =  actual  rainfall 

5  =:  maximum  retention  including  the  initial  abstraction. 

Infiltration,  infiltration  is  equal  to  the  difference  between  the  daily  precipitation, 
the  sum  of  the  change  in  the  surface  storage  of  precipitation  (snow),  the  daily  runoff, 
and  the  surface  evaporation.  Thus  the  net  daily  infiltration,  INj,  is  given  by: 

IN,  =  P,  -  Q,  -  ESS, 


where 

ESSi  =  surface  water  evaporation  on  day  /,  inches. 
Qj  =  daily  runoff,  inches. 

P,  =  net  rainfall,  given  by: 

p,  =  PRE,  -  SNQ -  SNOi 


where: 

PRE,  =  actual  precipitation  on  day  /,  in  inches. 

SNO,  =  amount  of  snowwater  at  end  of  day  /,  inches. 

Evapotranspiration.  The  potential  evapotranspiration  calculation  is  also  adapted 
from  Knisei  (Reference  52),  and  is  computed  using  a  modified  Penman  method 
developed  by  Ritchie  (Reference  53).  The  potential  evapotranspiration  on  day  /,  E^, 
is  given  by: 


e  1.28  A,  H, 
^  ■  (A,^  G)  25.4 


where: 

H,  >•  net  soier  radiation  on  day  /,  langleys. 

G  s  psychrometric  constant  which  equals  0.68 

A,  «  slope  of  saturation  vapor  pressure  curve  on  day  i,  computed  from: 

^  _  5304  _ ^21. 255-5304^ 

exp  - J 

where: 

T,  s  mean  temperature  in  on  day  /. 


The  daily  potential  evapotranspiration  demand  is  applied  first  on  water  available 
on  the  surface.  Any  demand  in  excess  of  the  surface  water  is  exerted  on  the  soil 
column  in  the  forms  of  soil  evaporation  and  plant  transpiration.  The  potential  soil 
evaporation,  ES^,  is  given  by: 


3-0.4LAI, 


where  L4/y  is  the  leaf  area  index  on  day  i,  or  the  winter  cover  factor  in  non>growing 
seasons.  Soil  evaporation  proceeds  at  this  rate  when  evaporation  is  not  limited  by 
transmission  of  water  to  the  surface.  Again  following  Knisei  (Reference  52),  this  limit 
is  given  by: 

M  9(a.-3)»« 

25.4 


where  a,  ~  soil  transmissivity  parameter  for  evaporation,  (mm/dayf  '^. 

After  reaching  this  limit,  soil  evaporation  proceeds  at  a  stage-two  rate,  £52,, 
given  by: 


ES^  . 


where  f,  =  days  since  stage  one  evaporation  ended. 

The  potential  plant  transpiration,  EP^,  is  computed  from: 

^  LAI, 

EP..  =  .  ' 
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Actual  plant  transpiration  depends  on  soil  moisture  and  plant  transpiration 
demand,  where  the  plant  transpiration  demand,  EPDj,  is  equal  to  the  potential  plant 
transpiration  except  when  limited  by  low  soil  moisture  or  when  the  daily  total  of  the 
surface  evaporation,  soil  evaporation  and  plant  transpiration  exceeds  the  daily 
potential  evapotranspiration.  The  actual  plant  transpiration,  £P,-  is  then  given  by: 


EP,  =  EPq 


1.20  -  (4  EPP) 


SM,  -  WP 
FC-WP 


in  which: 

SM,  =  depth  weighted  soil  moisture  on  day  /, 

WP  =  depth  weighted  wilting  point, 

PC  -  depth  weighted  field  capacity. 

Full  details  of  the  weighting  procedure  are  given  in  Schroeder  et  a/.  (Reference 


Soil  Moisture  Storage.  The  HELP  model  uses  a  daily  time  interval  to  evaluate 


the  components  of  the  water  balance  equation.  Soil  moisture  storage  is  then  given 
in  general  terms  as: 

SM,  = 

SM,.,  .  (  //V,  -  -  £7;  .  //V,.,  -  -  £7-,.,  ) 

where: 

SM,  = 

soil  moisture  storage  at  midday  / 

IN,  = 

infiltration  on  day  / 

PE,  = 

percolation  and  drainage  from  landfill  on  day  / 

ET,  = 

evapotranspiration  on  day  / 

In  application  of  the  model  the  vegetative  or  evaporative  zone  is  divided  into 
seven  segments.  Soil  water  is  then  distributed  among  these  segments  and  all 
underlying  layers,  with  the  equations  connected  by  vertical  drainage  terms.  The  model 
assumes  that  barrier  layers  always  remain  at  saturation.  After  distributing  the  water 
among  the  layers  the  model  checks  to  see  that  the  soi!  moisture  storage  does  not 
exceed  the  saturated  capacity.  Any  excess  above  this  amount  is  added  to  the  soil 
moisture  storage  of  the  layer  above,  or,  if  an  excess  in  the  topmost  segment,  added 
to  the  surface  runoff. 

Vertical  Flow.  The  model  assumes  that  the  soil  profile  consists  of  discrete 
segments  that  are  homogeneous  with  respect  to  the  hydraulic  conductivity,  the  total 
porosity  and  field  capacity.  The  rate  of  downward  flow  is  assumed  to  follow  Darcy's 
law.  However,  it  is  further  assumed  that  free  outflow  occurs  from  each  segment 
above  the  barrier  soil  layer,  in  which  case  the  rate  of  flow  equals  the  hydraulic 
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conductivity,  k.  This  assumption  is  valid  as  long  as  the  hydraulic  conductivities  of  the 
segments  above  the  barrier  soil  layer  are  similar,  or  increase  with  increasing  depths 
of  the  segments.  The  effective  hydraulic  conductivity,  k^,  is  a  function  of  the 
saturated  hydraulic  conductivity,  k„  and  the  soil  moisture  content,  defined  through: 


-K.[ 


SM,  -  MDC 
UL  -  MDC , 


where  MDC  is  the  minimum  soil  water  content  for  drainage  to  occur. 

The  routing  of  moisture  from  segment  to  segment  is  done  using  a  routing 
procedure  computed  at  the  mid-point  of  the  time  interval,  proceeding  sequentially  from 
the  top  segment  to  the  bottom  segment.  The  model  then  solves  simultaneously  for 
drainage  rate  and  soil  moisture.  If  the  moisture  content  of  a  segment  is  greater  than 
its  total  porosity  the  excess  is  backed  up  into  the  segment  above  it. 

After  convergence  is  obtained  for  the  segments  above  a  barrier  layer,  the 
hydraulic  head  may  be  computed  on  the  barrier.  Flow  through  the  barrier  layer  is  then 
also  computed  using  Darcy's  Law,  as: 

_  _  .  TH^TS(n.l) 

^  ~  •  TS(n-^l) 


where: 


k, 

TH  = 
TSfn+  1) 


saturated  hydraulic  conductivity  of  the  barrier  layer 
total  head  on  the  barrier  layer 
=  thickness  of  the  barrier  soil  layer. 


Lateral  Flow.  The  lateral  drainage  procedure  is  based  on  the  Boussinesq 
equation,  which  is  unsteady  and  non-lin'^  r: 

t  0h  _  IX  8f/t,„_\5hi.o 


where: 


dimensionless  drainable  porosity 
time  in  days 

gravitational  head  in  inches 

effective  saturated  lateral  hydraulic  conductivity,  inches/day 
lateral  position  in  the  direction  of  drainage 
dimensionless  slope 

recharge  flux  in  inches/day  perpendicular  to  the  direction  of  flow 


With  the  selection  of  a  small  time  step,  the  steady  state  assumption  is  made, 
d/i/df  =  0.  The  equation  is  then  linearized  following  the  form  given  by  Skaggs 
(Reference  54),  but  incorporating  an  additional  correction  factor  developed  by 
Schroeder  ef  aL  (Reference  50).  This  yields  the  approximate  relationship  for  lateral 
drainage,  QLAT,  as: 

QLAT  = 


where: 

L  <=  lateral  distance  from  the  crest  to  the  drain,  inches 

y  -  average  thickness  of  water  profile  above  barrier  soil  layer  between 

x=0  and  x=L,  inches. 

The  vertical  and  lateral  flow  routines  are  then  linked  under  two  assumptions: 
(1)  steady-state  conditions  hold  such  that  change  in  head  is  not  a  function  of  time, 
and  (2)  the  drainage  rate  estimated  at  the  mid-point  of  the  time  interval  is  effective 
throughout  the  time  interval.  These  assumptions  are  valid  only  if  the  computational 
time  interval  is  sufficiently  small  so  that  there  is  little  change  in  head.  The  model  is 
set  to  use  a  time  step  appropriate  for  most  common  conditions.  The  authors  contend 
that  four  equal  time  steps  per  day  yields  acceptable  accuracy  for  heads  less  than  30 
inches.  Finally,  the  model  must  perform  a  convergence  to  ensure  that  the  drainage 
rate  from  the  bottom  of  the  profile  is  equal  to  the  sum  of  lateral  flow  and  vertical 
percolation.  The  results  are  converged  to  the  5%  level  by  an  iterative  scheme 
commencing  with  an  a  priori  estimate  of  drainage  rate  from  the  bottom  profile' 
segment.  This  estimate  is  updated  until  convergence  within  the  5%  level  is  obtained. 


0.510  +  0.  00205  aL)  Kv 

L* 
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SECTION  III 


APPLICATIONS 


By  searching  the  Air  Force  IRPIMS  data  base,  it  was  determined  that  Hill  AFB 
(located  about  25  miles  north  of  Salt  Lake  City  and  5  miles  south  of  Ogden,  Utah)  had 
several  sites  that  are  good  candidates  for  mathematical  modeling.  This  was  confirmed 
after  a  2-day  visit  to  Hill  AFB  (June  25-26,  1990)  by  the  first  writer,  as  an  Air  Force- 
UES  Summer  Faculty  Fellow  at  Brooks  AFB,  Texas.  Hydrogeologic  and  contaminant 
transport  data  were  requested  for  these  sites:  a  limited  amount  of  data  was  obtained 
for  Operable  Unit  3  (0U3)  at  Hill  AFB  (References  55  and  56). 


A.  APPLICATION  OF  ODAST 

The  one-dimensional  analytical  model  ODAST  was  applied  by  James  M. 
Montgomery  (Reference  56)  to  illustrate  the  relationship  between  OU3  sources  and 
contaminant  migration.  The  model  assumes  that  the  aquifer  is  homogeneous  and 
isotropic  with  steady-state  uniform  groundwater  flow  at  constant  velocity,  it 
calculates  the  contaminant  concentration  at  any  time,  at  any  longitudinal  distance 
from  the  source  —  based  on  the  length  of  time  the  contaminant  was  injected  into  the 
groundwater,  and  it  adjusts  the  concentration  for  dispersion  and  retardation.  The 
resulting  concentration  profile  is  presented  in  Figure  6. 


B.  APPLICATION  OF  MOC  AND  REMEDIATION  OPTIMIZATION 

The  two-dimensional  flow  and  transport  model  developed  by  Konikow  and 
Bredehoeft  (Reference  21 )  was  used  to  simulate  groundwater  transport  around  OU3, 
Hill  AFB,  Utah.  This  model  has  been  updated  for  simulating  the  transport  of  non¬ 
conservative  contaminants.  The  area  of  the  regional  model  comprises  all  the  potential 
sources  in  OU3  and  a  substantial  area  downgradient:  with  uniformly  spaced  cells  at 
250-foot  intervals  (38  columns  and  39  rows).  All  recharges,  discharge,  and  leakage 
were  determined  for  each  cell  based  on  flow  conditions  established  by  the  regional 
groundwater  flow  model  MODFLOW.  Values  used  for  permeability  and  aquifer 
thickness  and  other  hydrogeologic  properties  were  the  same  as  those  for  the  regional 
groundwater  flow  model. 

Figure  7  depicts  an  isoconcentration  map  of  simulated  TCE  concentration,  with 
contour  intervals  from  1  ;/g/l  to  21  >7g/l,  for  a  retardation  factor  of  1.  Figure  8  is  a 
three-dimensional  view  of  the  simulated  TCE  concentration  contours. 
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Concentrafon  Profile  :  ODAST  model 
OU3.  Hill  AFB,  Utah 


Distance  From  ESE-9.  ft 


Figur*  6.  Concentration  Profile  as  Predicted  by  ODAST 
Optimization 

This  example  considers  only  pump  and  treat  activities  as  a  possible  remediation 
strategy.  In  addition,  the  hydraulic  conductivity,  retardation  by  sorption  and  porosity 
are  modeled  as  random  variables.  Figure  9  presents  a  map  of  the  simulated  contours 
of  TCE  contamination  and  pumping  well  placement  at  Hill  AFB.  Figure  10  presents 
the  results  of  evaluating  the  chance  constraint  with  only  1 00  iterations  and  maximum 
allowable  concentrations  ranging  from  0.1  mg/I  to  9.0  mg/I,  for  5  different 
remediation  strategies.  The  results  shown  here  define  the  cumulative  distribution  of 
the  contaminant  concentration  for  this  situation.  For  comparison.  Figure  10  also 
presents  the  cumulative  distribution  when  no  remediation  activities  are  implemented 
at  the  site. 

For  any  maximum  allowable  contaminant  concentration  (MCL)  a  tradeoff 
relationship  can  be  developed  between  the  cost  of  remediation  and  the  probability  that 
the  remediation  strategy  will  not  succeed  in  containing  contaminant  concentrations 
below  5  mg/I.  Figure  1 1  presents  this  tradeoff  relationship.  Relationships  similar  to 
those  presented  in  Figures  10  and  1 1  can  be  used  by  the  Air  Force  Technical  Project 
Managers  as  an  aid  in  evaluating  and  ultimately  deciding  on  a  long-term  remediation 
strategy  for  implementation. 
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rigur*  7.  TCE  Contours  for  Operable  Unit  3,  Hill  AFB 
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Tigxurtt  t.  TCE  Concentration  Surface,  Operable  Unit  3, 
Hill  AFB,  Utah 


78 


1  3  6  7  9  !t  13  16  17  19  21  23  25  27  29  31  33  36  37 

39 

37 

36 

33 

31 

29 

27 

25 

23 

21 

19 

17 

16 

13 

11 

9 

7 

6 

3 

1 

1  3  6  7  9  11  13  15  17  19  21  23  26  27  29  31  33  36  37 


11  I  I  I  I  I  I  I  I  I  I  M  I  I  I  ril  I  I  11  I  I  I  I  I  I  I  I  I  I  I  I 


:  Pollutant 
Sources 

:  Monitoring 
WeU 

:  Pumping 
Wells 


BuUding  220 
Building  511,5141 


be<W 

A  B^ian  Pondd 


» I  «  I  «  » »  I »  » I »  I  '  ■  1 1 1 1  I  i  1 1 1  i  I  I  I  I  I  I 


39 

37 


33 

31 

29 

27 

25 

23 

21 

19 

17 

16 

13 

11 

9 

7 

6 

3 

1 


Figur*  9.  Remediation  Well  Placement,  OU3,  Hill  AFB 
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Cumulative  Distribution  of  Contaminant  Concentration 

C0U3^  Hill  AFB^  5  Remediation  Strategies^ 


Maximum  Al lowab I e  Contaminant  Concentration,  MCL 
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Figur*  10.  Cuaulative  Distributions  of  Contaainant  Concentration 
For  Remediation  Schemes 
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Figur*  11.  Renedlation  Cost/Probability  of  Failxire  Tradeoff 
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C.  APPLICATION  OF  RESSQ 


RESSQ  was  used  to  simulate  2-dimensional  advective  transport  under  the 
injection,  extraction,  and  natural-gradient  conditions  of  the  tracer  experiments.  The 
field  site  Is  in  the  Moffett  Naval  Air  Station,  Mountain  View,  California  (Reference  57). 
The  RESSQ  model  was  used  to  estimate:  (1)  the  areal  extent  of  the  injection  fluid 
front  that  develops  around  the  injection  well  and  observation  wells,  (2)  the  fluid 
residence  times  from  the  injection  well  to  the  observation  wells,  and  (3)  the  degree 
of  recovery  of  the  injected  fluid  at  the  extraction  well.  A  sketch  of  the  the  well  fields 
is  presented  in  Figure  12,  for  fluid  injection  at  a  rate  of  0.5  liter/min  at  three  wells. 


Figure  12.  Map  of  Well  Fields  at  Moffett  Site 


extraction  rate  of  8  liters/min,  regional  ground-water  flow  of  300  m/yr,  a  porosity  of 
0.35,  and  an  aquifer  thickness  of  1.2  meters. 
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The  results  indicate  that  it  is  advantageous  to  use  the  southern  leg  for  the 
biorestoration  experiments.  The  reasons  are  the  following:  (1)  the  injected  fluid 
supplying  the  nutrients  becomes  less  dispersed,  and  hence  a  more  dense  microbial 
population  can  be  stimulated;  (2)  by  injecting  upgradient,  the  injected  tracers  and 
chlorinated  hydrocarbons  can  be  most  effectively  recovered  at  the  extraction  well. 
Figure  1 3  illustrates  the  streamlines  produced  by  a  graphics  post-processor  developed 
specifically  for  the  RESSQ  model. 


Figure  13.  RESSQ  Streamline  Plot 


SECTION  IV 


CONCLUSIONS 


The  objective  of  this  research  is  to  develop  a  groundwater  quality  and 
remediation  modeling  advisory  system  for  use  in  investigating  possible  strategies  for 
the  cleanup  of  contamination  from  hazardous  substances,  pollutants  and  contaminants 
at  Air  Force  sites.  In  addition,  the  use  of  optimization  methods  is  explored  for 
determining  optimal  remediation  for  implementation  at  a  specific  site.  The  modeling 
system  was  tested  successfully  for  Operable  Unit  3,  Hill  AFB,  Utah.  The  uncertainty 
associated  with  groundwater  flow  and  mass  transport  models  was  accounted  for,  not 
only  in  the  plume  prediction  itself,  but  also  explicitly  in  the  optimization  scheme.  The 
Advisory  System  also  provides  guidance  to  the  less  experienceu  users  in  proper  model 
selection.  This  report  constitutes  a  technical  summary  of  the  methodology  and  several 
of  the  component  models.  Additional  models  and  refinements  to  the  existing  codes 
are  planned  for  the  second  phase  of  the  research  project,  as  well  as  the  preparation 
of  a  users  manual. 

Mathematical  models  are  abstractions  of  the  real  physical  system  and/or 
biochemical  processes.  Over  the  past  several  decades  many  different  models  for 
contaminant  transport  in  porous  media,  under  varying  conditions  and  assumptions, 
have  been  proposed  and  tested.  These  range  from  very  simple  models  to  very 
complex  models.  All  of  these  models,  regardless  of  the  complexity  of  the  solution 
method,  require  certain  assumptions  regarding  the  nature  of  the  transport  processes, 
and  therefore  can  provide  only  an  approximation  of  the  actual  spread  of  contaminants 
from  a  given  site  and  the  associated  risks  from  human  exposure  to  contaminated 
groundwater.  This  situation  presents  a  familiar,  yet  difficult  problem  to  the  analyst 
and  the  decision-makers.  Sufficient  data  on  the  hydrogeology  are  rarely,  if  ever, 
available  to  apply  the  most  complex  codes.  The  analyst  must,  whether  explicitly  or 
implicitly,  choose  a  transport  model  based  on  a  trade-off  between  the  presumed 
greater  accuracy  of  complex  models  and  the  less  onerous  data  requirements  and  easier 
application  of  simpler  models:  an  important  justification  for  the  development  of 
assistance  tools  in  modeling,  despite  the  concerns  that  any  of  these  models  could  still 
be  incorrectly  applied  by  some  users.  It  is  ultimately  the  responsibility  of  the  analyst 
to  assess  the  results. 

Several  meetings  with  Air  Force  Center  for  Environmental  Excellence  (AFCEE) 
and  Air  Force  Engineering  and  Services  Center  (AFESC)  staff  and  representatives  of 
the  Installation  Restoration  Program  from  other  major  commands  were  held  to  solicit 
comment  from  potential  users.  Their  valuable  suggestions  are  being  incorporated  to 
make  the  system  more  useful  and  responsive  to  Air  Force  needs.  There  is 
considerable  debate  among  the  modeling  community  worldwide  regarding  validation 
of  groundwater  models,  since  they  embody  scientific  hypotheses  still  under 
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investigation.  The  models  assembled  in  this  Advisory  System  within  a  decision¬ 
making  framework  are  still  very  useful  in  developing  a  conceptual  model  of  behavior 
of  contaminants,  testing  which  remedial  strategies  are  the  most  reasonable,  improving 
data  collection  in  the  field  by  identifying  gaps  in  the  database  and  identifying  the  most 
sensitive  parameters  —  even  if  the  user  is  unable  to  fully  validate,  or  only  partially 
calibrate  the  results  to  field  data. 

In  addition  to  aiding  in  the  choosing  of  an  appropriate  mathematical  model  for 
a  specific  site,  the  Advisory  System  is  currently  being  modified  to  determine  efficient 
or  optimal  remediation  strategies.  The  optimization  routine  evaluates  tradeoffs 
between  the  long-term  cost  of  remediation  and  the  probability  the  remediation  strategy 
will  fail.  The  development  of  an  efficient,  effective  and  reliable  remediation  strategy 
requires  a  clear  understanding  of  the  site  characteristics  and  the  remediation  actions 
implemented,  in  addition,  the  optimal  remediation  strategy  must  consider  trade-offs 
between  the  remediation  cost  and  the  reliability  of  the  remediation  strategy.  By 
investigating  these  trade-offs,  the  decision  maker  can  more  accurately  assess 
remediation  needs,  feasible  remediation  strategies  and  remediation  strategy 
effectiveness.  Thus,  the  complete  system  incorporates  deterministic,  stochastic  and 
optimization  models  within  a  user-friendly  framework. 

Long-term  remediation  costs  depend  on  specific  remediation  considerations  and 
actions.  Examples  of  possible  remediation  strategies  include  pulse  pumping  and 
treatment,  and  continuous  pumping  and  treatment.  Potential  cost  savings  are  realized 
by  varying  the  long-term  remediation  action.  The  reliability  of  the  long-term 
remediation  strategy  represents  the  likelihood  that  contaminant  concentrations  within 
the  groundwater  exceeds  specified  maximums  and  are  modeled  as  constraints.  These 
two  conflicting  goals  or  objectives  are  weighed  against  one  another  using  a  chance 
constrained  optimization  model  in  which  the  physical  constraints  are  originally 
expressed  as  probabilistic  statements. 

Using  this  methodology,  optimal  groundwater  remediation  strategies  are 
determined  by  minimizing  the  long-term  and  short-term  costs  associated  with  the  site 
remediation.  In  addition,  the  optimal  remediation  strategies  are  conditioned  on  the 
probability  that  the  contaminant  concentration  at  any  time  does  not  exceed 
pre-specified  maxima.  The  actual  concentrations  at  any  specified  coordinates  are 
calculated  by  solving  the  governing  differential  equations  for  ground  water  contaminant 
flow  in  which  key  site  characteristics  are  expressed  as  random  variables.  The 
resulting  optimization  model  is  solved  using  a  second  moment  formulation  combined 
with  Monte  Carlo  simulation.  An  example  oJ  how  the  optimization  of  the  remediation 
process  works  has  been  demonstrated  for  Operable  Unit  3,  Hill  AFS,  Utah  for  five 
different  remedial  pumping  schemes. 
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